{
"cells": [
{
"cell_type": "markdown",
"id": "9cc3b252-0afc-404f-a7c2-706d0a7e3c89",
"metadata": {
"editable": true,
"id": "9cc3b252-0afc-404f-a7c2-706d0a7e3c89",
"slideshow": {
"slide_type": ""
},
"tags": [
"hide-input"
]
},
"source": [
"### Asking scientific questions of models - Exercises & Answers"
]
},
{
"cell_type": "markdown",
"id": "a734dcce-afbf-409a-8c38-ce4a535f4ea5",
"metadata": {
"id": "a734dcce-afbf-409a-8c38-ce4a535f4ea5"
},
"source": [
"The exercises here are designed to get you comfortable using models to make predictions and having them answer questions of interest, as opposed to relying on a suite of tests picked from a flowchart."
]
},
{
"cell_type": "markdown",
"id": "0cde8377-c11f-4474-95f2-7e6d353ecba1",
"metadata": {
"id": "0cde8377-c11f-4474-95f2-7e6d353ecba1"
},
"source": [
"## Traditional approaches from a model-based perspective\n",
"To get things clear, lets do some standard approaches such as t-tests and ANOVAs from the perspective of a linear model. We won't interpret the coefficients, we'll just get the model to tell us the answer directly and compare it to the traditional answer.\n",
"\n",
"### a. Imports\n",
"Import `pandas`, `pingouin`, `statsmodels.formula.api`, `seaborn`, and also `marginaleffects`."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "83f7ae21-b1f1-4d9b-93b2-720c03bbe75b",
"metadata": {
"editable": true,
"executionInfo": {
"elapsed": 11,
"status": "ok",
"timestamp": 1723926881418,
"user": {
"displayName": "Alex Jones",
"userId": "11094282981700434339"
},
"user_tz": -60
},
"id": "83f7ae21-b1f1-4d9b-93b2-720c03bbe75b",
"slideshow": {
"slide_type": ""
},
"tags": [
"hide-input"
]
},
"outputs": [],
"source": [
"# Your answer here\n",
"import pandas as pd\n",
"import pingouin as pg\n",
"import statsmodels.formula.api as smf\n",
"import seaborn as sns\n",
"import marginaleffects as me\n",
"\n",
"sns.set_style('whitegrid')"
]
},
{
"cell_type": "markdown",
"id": "9dc75c2e-25e4-4ed9-8f57-b60a682e038e",
"metadata": {
"id": "9dc75c2e-25e4-4ed9-8f57-b60a682e038e"
},
"source": [
"### b. Loading up data\n",
"We will continue our exploration of the 'Teaching Ratings' dataset here, and use `marginaleffects` to explore the consequences of our models.\n",
"\n",
"The data can be found here: https://vincentarelbundock.github.io/Rdatasets/csv/AER/TeachingRatings.csv\n",
"\n",
"Read it into a dataframe called `profs`, and show the top 5 rows."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "a87dbb65-80a5-4a29-a8b6-e26838ab5792",
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 224
},
"editable": true,
"executionInfo": {
"elapsed": 10,
"status": "ok",
"timestamp": 1723926881418,
"user": {
"displayName": "Alex Jones",
"userId": "11094282981700434339"
},
"user_tz": -60
},
"id": "a87dbb65-80a5-4a29-a8b6-e26838ab5792",
"outputId": "efe2f990-4fea-415e-f736-d6ae07fc0649",
"slideshow": {
"slide_type": ""
},
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
rownames
\n",
"
minority
\n",
"
age
\n",
"
gender
\n",
"
credits
\n",
"
beauty
\n",
"
eval
\n",
"
division
\n",
"
native
\n",
"
tenure
\n",
"
students
\n",
"
allstudents
\n",
"
prof
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
1
\n",
"
yes
\n",
"
36
\n",
"
female
\n",
"
more
\n",
"
0.289916
\n",
"
4.3
\n",
"
upper
\n",
"
yes
\n",
"
yes
\n",
"
24
\n",
"
43
\n",
"
1
\n",
"
\n",
"
\n",
"
1
\n",
"
2
\n",
"
no
\n",
"
59
\n",
"
male
\n",
"
more
\n",
"
-0.737732
\n",
"
4.5
\n",
"
upper
\n",
"
yes
\n",
"
yes
\n",
"
17
\n",
"
20
\n",
"
2
\n",
"
\n",
"
\n",
"
2
\n",
"
3
\n",
"
no
\n",
"
51
\n",
"
male
\n",
"
more
\n",
"
-0.571984
\n",
"
3.7
\n",
"
upper
\n",
"
yes
\n",
"
yes
\n",
"
55
\n",
"
55
\n",
"
3
\n",
"
\n",
"
\n",
"
3
\n",
"
4
\n",
"
no
\n",
"
40
\n",
"
female
\n",
"
more
\n",
"
-0.677963
\n",
"
4.3
\n",
"
upper
\n",
"
yes
\n",
"
yes
\n",
"
40
\n",
"
46
\n",
"
4
\n",
"
\n",
"
\n",
"
4
\n",
"
5
\n",
"
no
\n",
"
31
\n",
"
female
\n",
"
more
\n",
"
1.509794
\n",
"
4.4
\n",
"
upper
\n",
"
yes
\n",
"
yes
\n",
"
42
\n",
"
48
\n",
"
5
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" rownames minority age gender credits beauty eval division native \\\n",
"0 1 yes 36 female more 0.289916 4.3 upper yes \n",
"1 2 no 59 male more -0.737732 4.5 upper yes \n",
"2 3 no 51 male more -0.571984 3.7 upper yes \n",
"3 4 no 40 female more -0.677963 4.3 upper yes \n",
"4 5 no 31 female more 1.509794 4.4 upper yes \n",
"\n",
" tenure students allstudents prof \n",
"0 yes 24 43 1 \n",
"1 yes 17 20 2 \n",
"2 yes 55 55 3 \n",
"3 yes 40 46 4 \n",
"4 yes 42 48 5 "
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Your answer here\n",
"# Read in dataset\n",
"profs = pd.read_csv('https://vincentarelbundock.github.io/Rdatasets/csv/AER/TeachingRatings.csv')\n",
"profs.head()"
]
},
{
"cell_type": "markdown",
"id": "76eeeca9-97b3-463a-b485-5878d5b0124d",
"metadata": {
"id": "76eeeca9-97b3-463a-b485-5878d5b0124d"
},
"source": [
"### c. The t-test as a marginal effect\n",
"We will recreate a t-test with model-based predictions.\n",
"\n",
"First, conduct a t-test with `pingouin`, comparing the evaluation score of male and female professors."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "80728e27-f1d7-4b6d-8e87-9d9a7676463d",
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 114
},
"editable": true,
"executionInfo": {
"elapsed": 8,
"status": "ok",
"timestamp": 1723926881418,
"user": {
"displayName": "Alex Jones",
"userId": "11094282981700434339"
},
"user_tz": -60
},
"id": "80728e27-f1d7-4b6d-8e87-9d9a7676463d",
"outputId": "61160ca2-552d-4cbc-86a0-e110873cf56c",
"slideshow": {
"slide_type": ""
},
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
T
\n",
"
dof
\n",
"
alternative
\n",
"
p-val
\n",
"
CI95%
\n",
"
cohen-d
\n",
"
BF10
\n",
"
power
\n",
"
\n",
" \n",
" \n",
"
\n",
"
T-test
\n",
"
-3.266711
\n",
"
425.755804
\n",
"
two-sided
\n",
"
0.001176
\n",
"
[-0.27, -0.07]
\n",
"
0.305901
\n",
"
17.548
\n",
"
0.900288
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" T dof alternative p-val CI95% cohen-d \\\n",
"T-test -3.266711 425.755804 two-sided 0.001176 [-0.27, -0.07] 0.305901 \n",
"\n",
" BF10 power \n",
"T-test 17.548 0.900288 "
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Your answer here\n",
"# T-test with pingouin\n",
"pg.ttest(profs.query('gender == \"female\"')['eval'],\n",
" profs.query('gender == \"male\"')['eval']\n",
" )"
]
},
{
"cell_type": "markdown",
"id": "7af5a6c2-aa6b-4748-a1ed-194884866b78",
"metadata": {
"id": "7af5a6c2-aa6b-4748-a1ed-194884866b78"
},
"source": [
"Now fit a regression model with `statsmodels` predicting evaluations from gender. Call the model `ttest`. Check the summary, and remember the coefficient will equal the mean difference, which we can check our predictions against."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "69dadb4a-7a67-4d8a-8735-fad9245542c8",
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 253
},
"editable": true,
"executionInfo": {
"elapsed": 7,
"status": "ok",
"timestamp": 1723926881418,
"user": {
"displayName": "Alex Jones",
"userId": "11094282981700434339"
},
"user_tz": -60
},
"id": "69dadb4a-7a67-4d8a-8735-fad9245542c8",
"outputId": "dc7f5e9e-ef3a-46b3-dc75-a56ee5b415dc",
"slideshow": {
"slide_type": ""
},
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"
OLS Regression Results
\n",
"
\n",
"
Dep. Variable:
eval
R-squared:
0.022
\n",
"
\n",
"
\n",
"
Model:
OLS
Adj. R-squared:
0.020
\n",
"
\n",
"
\n",
"
No. Observations:
463
F-statistic:
10.56
\n",
"
\n",
"
\n",
"
Covariance Type:
nonrobust
Prob (F-statistic):
0.00124
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
coef
std err
t
P>|t|
[0.025
0.975]
\n",
"
\n",
"
\n",
"
Intercept
3.9010
0.039
99.187
0.000
3.824
3.978
\n",
"
\n",
"
\n",
"
gender[T.male]
0.1680
0.052
3.250
0.001
0.066
0.270
\n",
"
\n",
"
Notes: [1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
],
"text/latex": [
"\\begin{center}\n",
"\\begin{tabular}{lclc}\n",
"\\toprule\n",
"\\textbf{Dep. Variable:} & eval & \\textbf{ R-squared: } & 0.022 \\\\\n",
"\\textbf{Model:} & OLS & \\textbf{ Adj. R-squared: } & 0.020 \\\\\n",
"\\textbf{No. Observations:} & 463 & \\textbf{ F-statistic: } & 10.56 \\\\\n",
"\\textbf{Covariance Type:} & nonrobust & \\textbf{ Prob (F-statistic):} & 0.00124 \\\\\n",
"\\bottomrule\n",
"\\end{tabular}\n",
"\\begin{tabular}{lcccccc}\n",
" & \\textbf{coef} & \\textbf{std err} & \\textbf{t} & \\textbf{P$> |$t$|$} & \\textbf{[0.025} & \\textbf{0.975]} \\\\\n",
"\\midrule\n",
"\\textbf{Intercept} & 3.9010 & 0.039 & 99.187 & 0.000 & 3.824 & 3.978 \\\\\n",
"\\textbf{gender[T.male]} & 0.1680 & 0.052 & 3.250 & 0.001 & 0.066 & 0.270 \\\\\n",
"\\bottomrule\n",
"\\end{tabular}\n",
"%\\caption{OLS Regression Results}\n",
"\\end{center}\n",
"\n",
"Notes: \\newline\n",
" [1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
],
"text/plain": [
"\n",
"\"\"\"\n",
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: eval R-squared: 0.022\n",
"Model: OLS Adj. R-squared: 0.020\n",
"No. Observations: 463 F-statistic: 10.56\n",
"Covariance Type: nonrobust Prob (F-statistic): 0.00124\n",
"==================================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"----------------------------------------------------------------------------------\n",
"Intercept 3.9010 0.039 99.187 0.000 3.824 3.978\n",
"gender[T.male] 0.1680 0.052 3.250 0.001 0.066 0.270\n",
"==================================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
"\"\"\""
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Your answer here\n",
"# Fit the model\n",
"ttest = smf.ols('eval ~ gender', data=profs).fit()\n",
"ttest.summary(slim=True)"
]
},
{
"cell_type": "markdown",
"id": "309098eb-c5fa-4346-8302-2c8449632495",
"metadata": {
"id": "309098eb-c5fa-4346-8302-2c8449632495"
},
"source": [
"Next, use `marginaleffects` to create a datagrid that will give predictions for female and male professors, and pass it to `me.predictions` to make the predictions. Examine the values."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "e68cdbde-8412-462d-b45e-fcf709251814",
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 213
},
"editable": true,
"executionInfo": {
"elapsed": 370,
"status": "ok",
"timestamp": 1723926881781,
"user": {
"displayName": "Alex Jones",
"userId": "11094282981700434339"
},
"user_tz": -60
},
"id": "e68cdbde-8412-462d-b45e-fcf709251814",
"outputId": "5c953256-76f9-4720-c955-ef3df126a751",
"slideshow": {
"slide_type": ""
},
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"text/html": [
"
"
],
"text/plain": [
"shape: (1, 8)\n",
"┌───────────────┬──────────┬───────────┬──────┬─────────┬──────┬────────┬───────┐\n",
"│ Term ┆ Estimate ┆ Std.Error ┆ z ┆ P(>|z|) ┆ S ┆ 2.5% ┆ 97.5% │\n",
"│ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- │\n",
"│ str ┆ str ┆ str ┆ str ┆ str ┆ str ┆ str ┆ str │\n",
"╞═══════════════╪══════════╪═══════════╪══════╪═════════╪══════╪════════╪═══════╡\n",
"│ Row 1 - Row 2 ┆ 0.168 ┆ 0.0517 ┆ 3.25 ┆ 0.00115 ┆ 9.76 ┆ 0.0667 ┆ 0.269 │\n",
"└───────────────┴──────────┴───────────┴──────┴─────────┴──────┴────────┴───────┘\n",
"\n",
"Columns: term, estimate, std_error, statistic, p_value, s_value, conf_low, conf_high"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Your answer here\n",
"# Comparison is done via hypothesis\n",
"me.predictions(ttest, newdata=datagrid, hypothesis='pairwise')"
]
},
{
"cell_type": "markdown",
"id": "7af377ff-fbb7-4597-88ee-756f80eda804",
"metadata": {
"id": "7af377ff-fbb7-4597-88ee-756f80eda804"
},
"source": [
"### d. Carrying out an ANOVA with linear models and marginal effects\n",
"Lets now demonstrate how an ANOVA can be executed easily with a linear model and the examination of marginal effects.\n",
"\n",
"First, use `pinoguin` to carry out an ANOVA on teaching evaluations, using tenure and gender as the factors - that is, examine whether male and female professors differ in their evaluations depending on whether they have achieved tenure or not."
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "ece27b38-a2db-48f8-97c6-e8a636037065",
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 173
},
"editable": true,
"executionInfo": {
"elapsed": 8,
"status": "ok",
"timestamp": 1723926881781,
"user": {
"displayName": "Alex Jones",
"userId": "11094282981700434339"
},
"user_tz": -60
},
"id": "ece27b38-a2db-48f8-97c6-e8a636037065",
"outputId": "79d963c0-ca43-41c0-d431-4d41fcd2bd8b",
"slideshow": {
"slide_type": ""
},
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
Source
\n",
"
SS
\n",
"
DF
\n",
"
MS
\n",
"
F
\n",
"
p-unc
\n",
"
np2
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
gender
\n",
"
3.628914
\n",
"
1.0
\n",
"
3.628914
\n",
"
12.615338
\n",
"
0.000422
\n",
"
0.026749
\n",
"
\n",
"
\n",
"
1
\n",
"
tenure
\n",
"
2.829395
\n",
"
1.0
\n",
"
2.829395
\n",
"
9.835936
\n",
"
0.001821
\n",
"
0.020979
\n",
"
\n",
"
\n",
"
2
\n",
"
gender * tenure
\n",
"
4.187913
\n",
"
1.0
\n",
"
4.187913
\n",
"
14.558608
\n",
"
0.000154
\n",
"
0.030743
\n",
"
\n",
"
\n",
"
3
\n",
"
Residual
\n",
"
132.035435
\n",
"
459.0
\n",
"
0.287659
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" Source SS DF MS F p-unc np2\n",
"0 gender 3.628914 1.0 3.628914 12.615338 0.000422 0.026749\n",
"1 tenure 2.829395 1.0 2.829395 9.835936 0.001821 0.020979\n",
"2 gender * tenure 4.187913 1.0 4.187913 14.558608 0.000154 0.030743\n",
"3 Residual 132.035435 459.0 0.287659 NaN NaN NaN"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Your answer here\n",
"# A Pingouin ANOVA\n",
"pg.anova(data=profs, dv='eval', between=['gender', 'tenure'])"
]
},
{
"cell_type": "markdown",
"id": "b6b860ac-19c3-4e85-8a55-2613a24b9f2a",
"metadata": {
"id": "b6b860ac-19c3-4e85-8a55-2613a24b9f2a"
},
"source": [
"This suggests there is a main effect of gender, tenure and an interaction. Usually we'd need to do post-hoc tests to explore these. But we can rely on marginal effects for a simpler interpretation. First, fit a linear regression that is the same as the ANOVA, predicting evaluation measures from gender, tenure, and its interaction. Call it an `anova_model`."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "5f87872b-c485-41e2-818b-e52a298b17a2",
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 295
},
"editable": true,
"executionInfo": {
"elapsed": 7,
"status": "ok",
"timestamp": 1723926881782,
"user": {
"displayName": "Alex Jones",
"userId": "11094282981700434339"
},
"user_tz": -60
},
"id": "5f87872b-c485-41e2-818b-e52a298b17a2",
"outputId": "d6700cb9-298f-468c-b937-42fadb421a43",
"slideshow": {
"slide_type": ""
},
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"
OLS Regression Results
\n",
"
\n",
"
Dep. Variable:
eval
R-squared:
0.072
\n",
"
\n",
"
\n",
"
Model:
OLS
Adj. R-squared:
0.066
\n",
"
\n",
"
\n",
"
No. Observations:
463
F-statistic:
11.82
\n",
"
\n",
"
\n",
"
Covariance Type:
nonrobust
Prob (F-statistic):
1.80e-07
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
coef
std err
t
P>|t|
[0.025
0.975]
\n",
"
\n",
"
\n",
"
Intercept
3.8600
0.076
50.890
0.000
3.711
4.009
\n",
"
\n",
"
\n",
"
gender[T.male]
0.5362
0.106
5.047
0.000
0.327
0.745
\n",
"
\n",
"
\n",
"
tenure[T.yes]
0.0552
0.088
0.627
0.531
-0.118
0.228
\n",
"
\n",
"
\n",
"
gender[T.male]:tenure[T.yes]
-0.4610
0.121
-3.816
0.000
-0.699
-0.224
\n",
"
\n",
"
Notes: [1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
],
"text/latex": [
"\\begin{center}\n",
"\\begin{tabular}{lclc}\n",
"\\toprule\n",
"\\textbf{Dep. Variable:} & eval & \\textbf{ R-squared: } & 0.072 \\\\\n",
"\\textbf{Model:} & OLS & \\textbf{ Adj. R-squared: } & 0.066 \\\\\n",
"\\textbf{No. Observations:} & 463 & \\textbf{ F-statistic: } & 11.82 \\\\\n",
"\\textbf{Covariance Type:} & nonrobust & \\textbf{ Prob (F-statistic):} & 1.80e-07 \\\\\n",
"\\bottomrule\n",
"\\end{tabular}\n",
"\\begin{tabular}{lcccccc}\n",
" & \\textbf{coef} & \\textbf{std err} & \\textbf{t} & \\textbf{P$> |$t$|$} & \\textbf{[0.025} & \\textbf{0.975]} \\\\\n",
"\\midrule\n",
"\\textbf{Intercept} & 3.8600 & 0.076 & 50.890 & 0.000 & 3.711 & 4.009 \\\\\n",
"\\textbf{gender[T.male]} & 0.5362 & 0.106 & 5.047 & 0.000 & 0.327 & 0.745 \\\\\n",
"\\textbf{tenure[T.yes]} & 0.0552 & 0.088 & 0.627 & 0.531 & -0.118 & 0.228 \\\\\n",
"\\textbf{gender[T.male]:tenure[T.yes]} & -0.4610 & 0.121 & -3.816 & 0.000 & -0.699 & -0.224 \\\\\n",
"\\bottomrule\n",
"\\end{tabular}\n",
"%\\caption{OLS Regression Results}\n",
"\\end{center}\n",
"\n",
"Notes: \\newline\n",
" [1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
],
"text/plain": [
"\n",
"\"\"\"\n",
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: eval R-squared: 0.072\n",
"Model: OLS Adj. R-squared: 0.066\n",
"No. Observations: 463 F-statistic: 11.82\n",
"Covariance Type: nonrobust Prob (F-statistic): 1.80e-07\n",
"================================================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"------------------------------------------------------------------------------------------------\n",
"Intercept 3.8600 0.076 50.890 0.000 3.711 4.009\n",
"gender[T.male] 0.5362 0.106 5.047 0.000 0.327 0.745\n",
"tenure[T.yes] 0.0552 0.088 0.627 0.531 -0.118 0.228\n",
"gender[T.male]:tenure[T.yes] -0.4610 0.121 -3.816 0.000 -0.699 -0.224\n",
"================================================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
"\"\"\""
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Your answer here\n",
"# A linear model equivalent\n",
"anova_model = smf.ols('eval ~ gender * tenure', data=profs).fit()\n",
"anova_model.summary(slim=True)"
]
},
{
"cell_type": "markdown",
"id": "81931246-3836-4ad2-acc8-2a8001ca2971",
"metadata": {
"id": "81931246-3836-4ad2-acc8-2a8001ca2971"
},
"source": [
"With a fitted model, we can easily explore the implications via the predictions.\n",
"\n",
"First, make a datagrid that gives predictions for tenure and gender. Call it `anova_predmat`, and then use the model to predict those scores, storing them in a dataframe called `anova_predictions`."
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "cb771ac4-c3cc-4ac6-87f8-b6c46b38172c",
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 309
},
"editable": true,
"executionInfo": {
"elapsed": 7,
"status": "ok",
"timestamp": 1723926881782,
"user": {
"displayName": "Alex Jones",
"userId": "11094282981700434339"
},
"user_tz": -60
},
"id": "cb771ac4-c3cc-4ac6-87f8-b6c46b38172c",
"outputId": "b0babe63-dba5-42b7-83f7-8a0d60fb242c",
"slideshow": {
"slide_type": ""
},
"tags": [
"hide-input"
]
},
"outputs": [],
"source": [
"# Your answer here\n",
"# Prediction grid\n",
"anova_predmat = me.datagrid(anova_model,\n",
" tenure=['yes', 'no'],\n",
" gender=['male', 'female'])\n",
"\n",
"# Output\n",
"anova_predictions = me.predictions(anova_model, newdata=anova_predmat)"
]
},
{
"cell_type": "markdown",
"id": "dfbe945e-3969-4513-8cfe-3e9d9279abf8",
"metadata": {
"id": "dfbe945e-3969-4513-8cfe-3e9d9279abf8"
},
"source": [
"It is always sensible to plot predictions before we begin interpretin them. Use `seaborn` to create a line plot that illustrates the interaction. Any way you want is fine - as long as the estimate is on the y axis."
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "71c09c74-134f-4e35-a722-e9232800ce21",
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 159
},
"editable": true,
"executionInfo": {
"elapsed": 216,
"status": "ok",
"timestamp": 1723927069461,
"user": {
"displayName": "Alex Jones",
"userId": "11094282981700434339"
},
"user_tz": -60
},
"id": "os0rDvsqDUV-",
"outputId": "c7b8f145-dd1a-4cc4-8119-ec38128e7322",
"slideshow": {
"slide_type": ""
},
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjIAAAGsCAYAAADHSE33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABPxUlEQVR4nO3de1yUZf7/8dfMwHCQg6J4QFFDRaHMQ7hulnnWPFRrutVWVnaytNXMrEw7uZpmm7pZtm5p1v4693XLxMRMrXU9lxYpeD5iioICguMMM/fvj0kSj4MCNwPv5+NxP4R77sNnRuR+e93XfV0WwzAMRERERPyQ1ewCRERERC6VgoyIiIj4LQUZERER8VsKMiIiIuK3FGRERETEbynIiIiIiN9SkBERERG/FWB2AWXN4/FQWFiI1WrFYrGYXY6IiIj4wDAMPB4PAQEBWK3nb3ep9EGmsLCQ1NRUs8sQERGRS9CyZUvsdvt5X6/0QeZUimvZsiU2m63Ujut2u0lNTS3144qIiPiLsrwWnjr2hVpjoAoEmVO3k2w2W5kEjrI6roiIiL8oy2vhxbqFqLOviIiI+C0FGREREfFbCjIiIiLityp9HxlfeDwenE5nifZxu90AOBwO9ZHxgd1uv2iHLRERkZKq8kHG6XSya9cuPB5PifYzDIOAgAD27Nmj8Wl8YLVaueKKKy74CJ2IiEhJVekgYxgGv/76KzabjdjY2BK1GBiGwYkTJwgJCVGQuQiPx8OBAwf49ddfadiwoT4vEREpNVU6yBQWFlJQUEBMTAyhoaEl2vfUiIPBwcG6MPsgOjqaAwcOUFhYSGBgoNnliIhIJVGlOy2c6uei2x1l79RnfOozFxERKQ2mBhm3282gQYN45plnLrrt1q1badWqFWvWrCn1OtSiUvb0GYuISFkwNci88cYbrF+//qLbnThxglGjRuFwOMqhKhEREfEXpgWZVatWsXjxYnr27HnRbV966SW6d+9eDlWJiIiIPzGls29WVhZjx45l5syZzJ0794LbfvHFF+zZs4eJEycyc+bMSz7nufpmuN1uDMMoWkri1PYl3a8srVmzhnvvvZf09HSzSznLqc/Y7Xarn4yISCVx6vd5Wfxe9/WY5R5kPB4Po0ePZvDgwbRo0eKC2+7YsYNp06bx0UcfXfagc6mpqedcHxAQwIkTJ0o8jswpJ06cuJyyStXJkycBKCgoMLmSs508eRKXy1UhQ5aIiFwaj2Hg9pz/Glseyj3IzJo1C7vdzqBBgy643cmTJxk5ciTPPvssMTExl33ec00x7nA42LNnDyEhIQQHB5foeCUZR2bz5s1MmDCB9PR0GjZsyI033sinn37K0qVLWblyJdOmTWP37t3UqVOHIUOGcNNNNwHwzDPPEBQUxKFDh1i7di1RUVHcc8893HPPPQBkZmbywgsvsHbtWmrUqEGfPn0Aih4l37t3Ly+//DIbN24kNDSUm266iWHDhmG325k3bx4ffPABkZGRpKam8vzzzxedtyxYrVYCAwNp2rRpiT9rERGpWDweg/9sPMDfF2/F6XLx39GdCQ0u3SeA3W63TwGp3IPMl19+SWZmJklJSQBFHXiXLFlSrONvamoqu3fvZuzYsYwdO7Zo/SOPPMItt9zCiy++WKLznmuKcZvNhsViKVouxcX2PX78OA899BC333477733Hrt27eKRRx7BYrGwZcsWhg4dyquvvkq3bt346aefGDp0KDVq1KBjx45YLBbmzZvHrFmzePPNN/n8888ZP348N954I3Xq1OGJJ56gRo0afP/99+Tl5fHoo48W1VRQUMDgwYPp27cvr7/+OtnZ2QwfPhzDMBg1ahQWi4VNmzYxefJkZs2ahcfjKdMni059TmU51buIiJS91TuzmJC8mV8ycgFoFBmALSDAtN/t5R5kFi1aVOz7U49eT548udj6pKQkfv7552Lrmjdvzj//+U/at29ftkWWoqVLl2Kz2fjrX/+K1WqlefPmPPjgg8yePZuPP/6Ybt26FXV4btu2LbfddhsffPABHTt2BKB9+/Zcd911AAwYMIAXXniBvXv3UlhYyPr160lJSSEsLIywsDAee+wxhg0bBsDy5ctxOp088cQTWCwW6tWrx4gRIxg+fDijRo0CIDAwkFtuuUVzIImIyEXtPpLPpK/TSNl0CIDwoACGdmlC69BjBAWYdx2pcCP7tmnThpdeeombb77Z7FJKxcGDB4mJiSkWFmJjYwHIyMhg9erVRa1T4G1Ka9iwYdH30dHRRV+fGhHX4/Fw6JD3B+n0226n75eRkUF2djbt2rUrWmcYBi6Xi6ysrKJjK8SIiMiF5BS4mLF0G++t2o3LbWC1wJ3tGzKyezzVQwLYuHGjqfWZHmTObInZsGHDebfdsmVLWZdT6mJiYjhw4ACGYRTdujlw4AAAdevWpX///owfP75o+8zMTJ+ehKpbty4A+/bto0mTJoA3NJ3+esOGDYu1gB0/fpysrCyioqIADVInIiLn53J7+HDNXqYv2crRAhcAnZtH82yfBOLrhAMVY7R2/Xe8jHXt2hXDMPjnP/+J0+lk586dzJ49G4CBAweyYMECVqxYgcfjYffu3dx9993MmTPnoseNiYnh+uuvZ9KkSeTk5HD48GHeeOONote7dOlCfn4+77zzDk6nk9zcXJ5++mlGjhypACMiIudlGAZL0w9x4/TveWH+Jo4WuIivE8Z79/+BuYP/UBRiKgoFmTIWGhrKzJkz+fbbb/nDH/7AE088wXXXXUdgYCCtWrVi6tSpTJ06lXbt2nH33XfTtWvXoj4sF/Paa68RHh5Oly5dGDBgAB06dCh6LSwsjLlz57JmzRpuuOEGunfvjtVq5a233iqrtyoiIn4u/WAug2av5f6569lxOJ+a1exM7H8VC4d3pFN89MUPYAKLUZFGdCsDbrebjRs30rp163M+fr1r1y6uuOKKS3r8uqCggNDQ0Au2cBw9epSdO3dyzTXXFK3797//TXJyMh9//HHJ3owfu5zPWkREylZmnoNp32zlk3X78Bhgt1kZfH1jhnVpSkRw4Hn3u9A19nL5emzT+8hUdm63m3vvvZc333yTTp06sX//fj788EP+9Kc/mV2aiIhUcQ6Xm9krdjFz2Xbynd7+Ln2vrsczN7YgNirU5Op8oyBTxmrVqsX06dP5+9//zuOPP05ERAT9+/fngQceMLs0ERGpogzDYP5PB5iyaAsZx7wj1LdqEMlz/RJJahxlcnUloyBTDrp3765JL0VEpEL4Yc9RJiRvZsPeYwDERAbzdO8W3HR1DFar/z0MoiAjIiJSBew/WsAri7bw1U/eIUBC7TaGdm7CA9fHEWL33xHXFWREREQqsTyHi7eW7+CdFbtwFnqwWOC2a2IZ1TOe2hH+//CFgoyIiEgl5PYYfLJuH1O/2cKR404Aro2rybh+CVwZE2lydaVHQUZERKSS+e+2w0xMTiP9YB4AcbWq8WyfBLol1K50g6IqyIiIiFQS2zPzmJicxrIthwGIDAnk8e7NuKt9I+wmTuxYlhRkRERE/Fx2vpPpS7bywZq9uD0GAVYL91zbmOHdmlI91G52eWVKQaaK69q1K4899hi33nqr2aWIiEgJnSx08/7KPby+dBt5jkIAeiTWYUzvFsRFh5lcXflQkDmDYRiccF18Nk/DMChwuiGg8LLvN4YE2irdPUsRESk7hmGQsukgk75OZ09WAQCJ9SIY1y+BDk1qmVxd+VKQOY1hGAz85yp+2HO0XM+b1KgGnz1yrc9hZv/+/XTr1o1XXnmFf/zjHxw9epTevXszYMAAxo8fz759+7j66quZNm0adrudyZMns3btWjIzMwkPD+euu+7ikUceOeu4TqeTt956i/nz55OXl0erVq0YN24cjRo1Ku23LCIilyh1fw5/S97M2l3ZAESHBzG6V3MGtG2AzQ8HtLtcCjJn8Kcfge+++46FCxeyb98+/vSnP7F582befvttAgMDueOOO/jwww85cuQI+/fv5/PPPyc8PJzFixczfPhwevfufVZAmTZtGqtXr2bu3LnUrl2bt99+m/vvv5+FCxcSFBRk0rsUERGAgzkOpqSkM+/HDACCA6083DGOIZ2aUC2o6l7Oq+47PweLxcJnj1zr+62lghOEhoaYdmvp/vvvJyQkhPj4eKKjo+nfvz916tQBoHXr1mRkZPDkk09is9kICwvj4MGDRYEkMzOzWJAxDIOPP/6Y119/ndjYWACGDRvGp59+yvLly+nVq9dlvUcREbk0Bc5CZn23k1nf78Dh8gDQv019RvdqTkz1EJOrM5+CzBksFguh9ot/LIZhQKGNUHuAaf1bqlevXvS1zWYjIiKi6Hur1YphGGRlZTFx4kQ2b95MgwYNuOqqqwDweDzFjpWdnU1BQQEjRozAav39ET2Xy0VGRkbZvhERETmLx2Mwb0MGr6akcyj3JADtGtdgXN9EWsVWN7e4CkRBxo/5EqBGjBhB165dmT17NgEBARw9epRPP/30rO1q1KhBUFAQc+bMoXXr1kXrd+7cWdTKIyIi5WP1ziwmJG/ml4xcAGKjQhjTO4HeV9XVwyFnqJyj40iRvLw8goODsdlsZGdnM2HCBMDb0nI6q9XKwIEDee211zh48CAej4f//Oc/9OvXjz179phRuohIlbP7SD5D/r2eO/61ml8ycgkPCmBM7xZ8M7ITfVrWU4g5B7XIVHKTJk3i5ZdfZs6cOURGRtKnTx8SExPZunUr119/fbFtn376aWbMmMGdd97JsWPHiI2N5fXXXycxMdGk6kVEqoacAhczlm7jvVW7cbkNrBa4s31DRnaPp2aYHra4EIthGIbZRZQlt9vNxo0bad26NTZb8WnKHQ4Hu3bt4oorriA4uGQzgHo7+xYQGhqqhOyDy/msRUQqK5fbw4dr9jJ9yVaOFnhbyjs3j+bZPgnE1wk3ubqLu9A1tryOrRYZERGRcmYYBsu2ZDIxOY0dh/MBiK8Txti+iXSKjza5Ov+iICMiIlKO0g/mMmFBGiu2HwGgZjU7I3vEc0e7WAJs6rpaUgoyIiIi5SAzz8G0b7byybp9eAyw26wMvr4xw7o0JSI40Ozy/JaCjIiISBlyuNzMXrGLmcu2k+/0Drja9+p6PHNjC2KjQk2uzv8pyIiIiJQBwzCY/9MBpizaQsaxEwC0ahDJc/0SSWocZXJ1lYeCjIiISCn7Yc9RJiRvZsPeYwDERAbzdO8W3HR1DNYqOLFjWVKQERERKSX7jxYw+et0Fvz8KwChdhtDOzfhgevjCLGX7uPJ4qUgIyIicpnyHC7eWr6Dd1bswlnowWKB266JZVTPeGpHaOyssqQgIyIiconcHoNP1u1j6jdbOHLcCcC1cTUZ1y+BK2MiTa6ualCQERERuQT/3XaYiclppB/MAyCuVjWe7ZNAt4TaGvG9HGnkHT+0cOFCrr32Wq655hqWLVtWLufcv38/zZs3Z//+/eVyPhGRimp7Zh6D313LoNlrST+YR2RIIC/clMiix2+ge2IdhZhyphaZMxkGuAp82855AgIMuNwf2sDQEh3js88+o2/fvowbN+7yzisiIj7LzncyfclWPlizF7fHIMBq4Z5rGzO8W1Oqh9rNLq/KUpA5nWHAnF6wb81FN7UA1UrrvLF/hPsX+RRmBg4cyKZNm1i3bh3Lly9nzpw5vPzyy2zYsIHQ0FBuvvlmhg0bht1uZ968eXz++ee0atWK//u//8NqtTJs2DCCgoJ46623yM3NpW/fvowfPx6AHTt2MGXKFLZs2UJ2djYNGjRg9OjRdOnS5aw6jhw5wuTJk1m1ahUWi4WuXbvy1FNPERYWVlqfiohIhXCy0M37K/fw+tJt5DkKAeiRWIcxvVsQF63feWbTraWzVOwmwc8//5ykpCSGDBnC/Pnzue+++2jWrBnff/89H374IStXrmTGjBlF2//www/UqVOH1atXM3z4cCZNmsSaNWtYuHAhc+fO5fPPP2fdunUA/PWvfyU+Pp5vvvmG9evXc/311/Piiy+eVYPH42Ho0KFYrVZSUlL46quvyMzM5Pnnny+vj0FEpMwZhsGiX36l57TvmbgwjTxHIYn1Ivjwwfa8fU+SQkwFoRaZ01ks3pYRH24tGYZBQcEJQkNDLv9+aAlvLZ2yfPlynE4nTzzxBBaLhXr16jFixAiGDx/OqFGjAAgNDeXee+/FYrFw/fXX43a7eeCBBwgJCaFly5bUrl2bjIwM2rVrx6xZs6hTpw6GYZCRkUFERASHDh0667y//PILmzZt4t1336VaNW+71NNPP82NN97Ic889R40aNS7v8xARMVnq/hz+lryZtbuyAYgOD2J0r+YMaNsAmwa0q1AUZM5ksYDdh5tGhgGFFrBfWggpDRkZGWRnZ9OuXbvTyjJwuVxkZWUBUL169aKgZbV6G+AiIiKKtrdarXg8HgDS09MZOnQohw8fpkmTJkRFRWEYxlnn3b9/P263m06dOhVbb7fb2bdvn4KMiPitgzkOpqSkM+/HDACCA6083DGOIZ2aUC1Il8yKSH8rfqxu3bo0bNiQRYsWFa07fvw4WVlZREV55/HwtbXo0KFDjBgxgjfeeIOuXbsCkJKSwuLFi8953uDgYNasWYPN5h2p0ul0sm/fPho1anS5b0tEpNwVOAuZ9d1OZn2/A4fL+5+7/m3qM7pXc2Kqh5hcnVyI+sj4sS5dupCfn88777yD0+kkNzeXp59+mpEjR5b4dld+fj5ut5uQEO8/2O3bt/Pmm28C3pByuquvvppGjRoxefJk8vPzcTgcvPzyy9x333243e7SeXMiIuXA4zH4/If9dPn7cv7x7TYcLg/tGtfgy2HXMe321goxfkBBxo+FhYUxd+5c1qxZww033ED37t2xWq289dZbJT5WXFwcTz31FKNHj+aaa65hxIgRDBgwgMDAQLZu3Vps24CAAGbNmsWRI0fo2bMn119/PXv37uXdd98lKCiotN6eiEiZWr0zi5vfXMGTn/3EodyTxEaFMPOutnw65FpaxVY3uzzxkcU4VyeISsTtdrNx40Zat25ddBvkFIfDwa5du7jiiisIDi7ZXBjezr4FhIaGavAjH1zOZy0iUpp2H8ln0tdppGzyPswQHhTAY12bcm+HxgQHamLHkrjQNba8jq0+MiIiUiXkFLiYsXQb763ajcttYLXAne0bMrJ7PDXD1JrsrxRkRESkUnO5PXy4Zi/Tl2zlaIELgM7No3m2TwLxdcJNrk4ul4KMiIhUSoZhsGxLJhOT09hxOB+A+DphjO2bSKf4aJOrk9KiICMiIpVO+sFcJixIY8X2IwDUrGZnZI947mgXS4BNz7lUJgoycM5B36R06TMWkfKQmedg2jdb+WTdPjwG2G1WBl/fmGFdmhIRHGh2eVIGqnSQOX0wt1Pjp0jZODUWTWn3ahcRAXC43MxesYuZy7aT7/SOZ9X36no8c2MLYqNCTa5OylKVDjIBAQGEhoZy+PBhAgMDi4bw94VhGJw8eRKr1arHry/C4/Fw+PBhQkNDCQio0j9yIlLKDMNg/k8HmLJoCxnHTgDQqkEkz/VLJKlxlMnVSXmo0leVUxMt7tq1iz179pRo31NzGgUGBirI+MBqtdKwYUN9ViJSan7Yc5QJyZvZsPcYADGRwTzduwU3XR2DVRM7VhlVOsiAd6LDZs2anTUM/8W43W7S09Np2rSpbpf4wG63l6jFS0TkfPYfLWDy1+ks+PlXAELtNh7t1IQHO8YRYtfv46qmygcZ8LYWlHS02VNzCgUHByvIiIiUgzyHi7eW7+CdFbtwFnqwWOC2a2IZ1TOe2hEaMbyqUpAREZEKze0x+GTdPqZ+s4Ujx72t59fG1WRcvwSujIk0uToxm6lBxu12c99991G/fn0mT5581usej4c333yTzz//nNzcXBo0aMCjjz5Knz59TKhWRETK23+3HWZichrpB/MAiKtVjWf7JNAtobb63AlgcpB54403WL9+PfXr1z/n6x988AFffPEF//73v2nYsCHLli1j6NChXHXVVTRs2LCcqxURkfKyPTOPiclpLNtyGIDIkEAe796Mu9o3wh6g/nbyO9OCzKpVq1i8eDE9e/Y87zZ33XUXAwYMIDQ0FKfTSXZ2NiEhIZc0e/KpPi2l5dTxSvu4IiJVWXa+k9e/3c6H6/bh9hgEWC0M+mNDHuvShOqhdsDQ790KpCyvhb4e05Qgk5WVxdixY5k5cyZz584973ZWq5XQ0FBWrFjBQw89hGEYjBkzhtq1a5f4nKmpqZdRcfkfV0SkKnG5Db7eXsBnaccpcHlHAm8XE8Q9V4cTE36S3Vs3m1yhXIiZ18JyDzIej4fRo0czePBgWrRo4dM+f/jDH0hNTWXdunUMHTqU6OjoEveTadmyZak+XeR2u0lNTS3144qIVCWGYbB48yEmf7uVvdkFACTWC2dM7xZ0aFLT5OrkYsryWnjq2BdT7kFm1qxZ2O12Bg0a5PM+drsdgGuvvZZbbrmFr776qsRBxmazlUngKKvjiohUdqn7c/hb8mbW7soGIDo8iNG9mjOgbQNsGtDOr5h5LSz3IPPll1+SmZlJUlISAA6HA4AlS5awfv36YtueepLpmWeeKVrndDqpXr16+RQrIiKl7mCOgykp6cz7MQOA4EArD3eMY0inJlQL0qggUjLl/hOzaNGiYt+fCinnevw6KSmJJ598km7dunHNNdewfPlyFi5cyJw5c8qlVhERKT0FzkJmfbeTWd/vwOHyANC/TX1G92pOTHVN3CuXpsJF3zZt2vDSSy9x88030717d8aNG8e4ceM4cuQIjRs3ZsaMGbRt29bsMkVExEcej8G8DRm8mpLOodyTACQ1qsFz/RJpFVvd3OLE75keZM5sidmwYUOx7wcOHMjAgQPLsyQRESklq3dmMSF5M79k5AIQGxXCmN4J9L6qrga0k1JhepAREZHKZ/eRfCZ9nUbKpkMAhAcF8FjXptzboTHBgXpAQkqPgoyIiJSanAIXM5Zu471Vu3G5DawWuLN9Q0Z2j6dmWJDZ5UklpCAjIiKXzeX28OGavUxfspWjBS4AOjeP5tk+CcTXCTe5OqnMFGREROSSGYbB0vRMJi5MY+fhfADi64Qxtm8ineKjTa5OqgIFGRERuSRpv+YyMTmNFduPAFCzmp2RPeK5o10sATZN7CjlQ0FGRERKJDPPwbRvtvLJun14DLDbrAy+vjHDujQlIjjQ7PKkilGQERERnzhcbmav2MXMZdvJd3pnJu57dT2eubEFsVGhJlcnVZWCjIiIXJBhGMz/6QBTFm0h49gJAFo1iOS5fokkNY4yuTqp6hRkRETkvH7Yc5QJyZvZsPcYADGRwTzduwU3XR2DVRM7SgWgICMiImfZf7SAyV+ns+DnXwEItdt4tFMTHuwYR4hdA9pJxaEgIyIiRfIcLt5avoN3VuzCWejBYoHbrollVM94akcEm12eyFkUZEREBLfH4JN1+5j6zRaOHHcCcG1cTcb1S+DKmEiTqxM5PwUZEZEq7r/bDjMxOY30g3kAxNWqxrN9EuiWUFsTO0qFpyAjIlJFbc/MY2JyGsu2HAYgMiSQx7s34672jbAHaEA78Q8KMiIiVUx2vpPpS7bywZq9uD0GAVYL91zbmOHdmlI91G52eSIloiAjIlJFnCx08/7KPby+dBt5jkIAeiTWYUzvFsRFh5lcncilUZAREankDMMgZdNBJn2dzp6sAgAS60Uwrm8CHZrWMrk6kcujICMiUoml7s/hb8mbWbsrG4Do8CBG92rOgLYNsGlAO6kEFGRERCqhgzkOpqSkM+/HDACCA6083DGOIZ2aUC1Iv/ql8tBPs4hIJVLgLGTWdzuZ9f0OHC4PAP3b1Gd0r+bEVA8xuTqR0qcgIyJSCXg8BvM2ZPBqSjqHck8CkNSoBs/1S6RVbHVzixMpQwoyIiJ+bvXOLCYkb+aXjFwAYqNCGNM7gd5X1dWAdlLpKciIiPip3UfymfR1GimbDgEQHhTAY12bcm+HxgQHamJHqRoUZERE/ExOgYsZS7fx3qrduNwGVgvc2b4hI7vHUzMsyOzyRMqVgoyIiJ9wuT18uGYv05ds5WiBC4DOzaN5tk8C8XXCTa5OxBwKMiIiFZxhGCxNz2TiwjR2Hs4HoFntMMb1S6RTfLTJ1YmYS0FGRKQCS/s1l4nJaazYfgSAmtXsjOwRzx3tYgmwaWJHEQUZEZEKKDPPwbRvtvLJun14DLDbrAy+vjHDujQlIjjQ7PJEKgwFGRGRCsThcjN7xS5mLttOvtMNQN+r6/HMjS2IjQo1uTqRikdBRkSkAjAMg/k/HWDKoi1kHDsBQKsGkTzXL5GkxlEmVydScSnIiIiY7Ic9R5mQvJkNe48BEBMZzNO9W3DT1TFYNbGjyAUpyIiImGRfdgGvLEpnwc+/AhBqt/FopyY82DGOELsGtBPxhYKMiEg5y3O4mLl8B7NX7MJZ6MFigduuiWVUz3hqRwSbXZ6IX1GQEREpJ26PwSfr9jH1my0cOe4E4Nq4mozrl8CVMZEmVyfinxRkRETKwX+3HWZichrpB/MAiKtVjWf7JNAtobYmdhS5DAoyIiJlaHtmHhOT01i25TAAkSGBPN69GXe1b4Q9QAPaiVwuBRkRkTKQne9k+pKtfLBmL26PQYDVwj3XNmZ4t6ZUD7WbXZ5IpaEgIyJSik4Wunl/5R5eX7qNPEchAD0S6zCmdwviosNMrk6k8lGQEREpBYZhkLLpIJO+TmdPVgEAifUiGNc3gQ5Na5lcnUjlpSAjInKZUvfn8LfkzazdlQ1AdHgQo3s1Z0DbBtg0oJ1ImVKQERG5RAdzHExJSWfejxkABAdaebhjHEM6NaFakH69ipQH/UsTESmhAmchs77byazvd+BweQDo36Y+o3s1J6Z6iMnViVQtCjIiIj7yeAzmbcjg1ZR0DuWeBCCpUQ3G9UukdWx1c4sTqaIUZEREfLB6ZxYTkjfzS0YuALFRIYzpnUDvq+pqQDsREynIiIhcwO4j+Uz6Oo2UTYcACA8K4LGuTbm3Q2OCAzWxo4jZFGRERM4hp8DFjKXbeG/VblxuA6sF7mzfkJHd46kZFmR2eSLyGwUZEZHTuNwePlyzl+lLtnK0wAVA5+bRPNsngfg64SZXJyJnUpAREcE7oN3S9EwmLkxj5+F8AJrVDmNcv0Q6xUebXJ2InI+CjIhUeWm/5jIxOY0V248AULOanZE94rmjXSwBNk3sKFKRKciISJWVmedg6uKtfLp+Hx4D7DYrg69vzLAuTYkIDjS7PBHxgYKMiFQ5Dpeb2St2MXPZdvKdbgD6Xl2PZ25sQWxUqMnViUhJKMiISJVhGAbzfzrAlEVbyDh2AoBWDSJ5rl8iSY2jTK5ORC6FqUHG7XZz3333Ub9+fSZPnnzObT766CPmzp1LZmYmtWvX5p577uGuu+4q50pFxN/9sOcoE5I3s2HvMQDqRQbz9I0tuLlVDFZN7Cjit0wNMm+88Qbr16+nfv3653x9yZIlTJ06lbfffptWrVqxceNGHn74YWrVqkWvXr3KuVoR8Uf7sgt4ZVE6C37+FYBQu41HOzXhwY5xhNg1oJ2IvzMtyKxatYrFixfTs2fP825z6NAhHnroIVq3bg1AmzZtaN++PevWrVOQEZELynO4mLl8B7NX7MJZ6MFigduuiWVUz3hqRwSbXZ6IlBJTgkxWVhZjx45l5syZzJ0797zbnXkLKSsri3Xr1jFmzJgSn9Ptdpd4H1+OV9rHFZHL4/YYfLp+P9OWbCMr3wnAH+OiGNu7BYkxEd5t9O9WpFSU5bXQ12OWe5DxeDyMHj2awYMH06JFC5/3O3z4MEOGDOGqq66iX79+JT5vampqifcx87giUnI/HTrJ3J/y2JtTCEBMmI17WoWTVC8QZ+ZONmaaXKBIJWXmtbDcg8ysWbOw2+0MGjTI5302btzIiBEjSEpKYtKkSQQElLzsli1bYrOV3v1wt9tNampqqR9XREpue+ZxJn29heVbjwIQGRLI8K5NuPMPDbEHaEA7kbJSltfCU8e+mHIPMl9++SWZmZkkJSUB4HA4AG/H3vXr15+1/eeff86ECRMYPnw4999//yWf12azlUngKKvjisjFZec7mb5kKx+s2YvbYxBgtXDPtY0Z3q0p1UPtZpcnUmWYeS0s9yCzaNGiYt8/88wzAOd8/DolJYUXX3yRt956i44dO5ZLfSJS8Z0sdPP+yj28vnQbeQ7vbaQeiXUY07sFcdFhJlcnIuWpwg2I16ZNG1566SVuvvlm3njjDdxuN8OHDy+2zU033cT48eNNqlBEzGIYBimbDjLp63T2ZBUAkFgvgnF9E+jQtJbJ1YmIGUwPMme2xGzYsKHo66+++qq8yxGRCip1fw5/S97M2l3ZAESHBzG6V3MGtG2ATQPaiVRZpgcZEZELOZjjYEpKOvN+zAAgONDKwx3jGNKpCdWC9CtMpKrTbwERqZAKnIXM+m4ns77fgcPlAaB/m/qM7tWcmOohJlcnIhWFgoyIVCgej8G8DRm8mpLOodyTACQ1qsG4fom0jq1ubnEiUuEoyIhIhbF6ZxYTkjfzS0YuALFRIYzpnUDvq+pisagfjIicTUFGREy3+0g+k75OI2XTIQDCgwJ4rGtT7u3QmOBAjdMkIuenICMipskpcDFj6TbeW7Ubl9vAaoE72zdkZPd4aoYFmV2eiPgBBRkRKXcut4cP1+xl+pKtHC1wAdApPpqxfROIrxNucnUi4k8UZESk3BiGwdL0TCYuTGPn4XwAmtUOY2zfBDo3r21ydSLijxRkRKRcpP2ay8TkNFZsPwJAzWp2RvaI5452sQTYNLGjiFwaBRkRKVOZeQ6mLt7Kp+v34THAbrMy+PrGDOvSlIjgQLPLExE/pyAjImXC4XIze8UuZi7bTr7TDUDfq+vxzI0tiI0KNbk6EaksFGREpFQZhsH8nw4wZdEWMo6dAKBVg0ie65dIUuMok6sTkcpGQUZESs0Pe44yIXkzG/YeA6BeZDBP39iCm1vFYNXEjiJSBhRkROSy7csu4JVF6Sz4+VcAQu02Hu3UhAc7xhFi14B2IlJ2FGRE5JLlOVzMXL6D2St24Sz0YLHAbdfEMqpnPLUjgs0uT0SqAAUZESkxt8fgk3X7mPrNFo4cdwJwbVxNxvVL4MqYSJOrE5Gq5JKCzP/+9z/+/e9/k5mZyaxZs5gzZw6jRo0iIEC5SKSy+++2w0xMTiP9YB4AcbWq8WyfBLol1NbEjiJS7kqcPL766ismTZrEn//8Z9atWwfA0qVLsVgsPPXUU6VeoIhUDNsz85iYnMayLYcBiAwJZES3Ztz9x0bYAzSgnYiYo8S/ff71r38xc+ZMRo4cidVqJTo6mlmzZrFgwYKyqE9ETJad7+T5L3+h1/T/smzLYQKsFu6/7gq+G92Z+6+/QiFGRExV4haZgwcP0qpVK4CiZuRGjRpRUFBQupWJiKlOFrp5f+UeXl+6jTxHIQA9EuswpncL4qLDTK5ORMSrxEGmcePGfPvtt3Tv3r1o3cqVK2nUqFGpFiYi5jAMg5RNB5n0dTp7srz/QUmsF8G4vgl0aFrL5OpERIorcZAZOXIkQ4cOpVu3bpw8eZIXX3yRBQsW8Nprr5VFfSJSjlL35/C35M2s3ZUNQHR4EKN7NWdA2wbYNKCdiFRAJQ4yHTp04OOPP+aTTz6hffv2eDwe5syZw9VXX10W9YlIOTiY42BKSjrzfswAIDjQysMd4xjSqQnVgvQ0oohUXCX+DTV79mweeOABXnjhhWLrp0+fzuOPP15adYlIOShwFjLru53M+n4HDpcHgP5t6jO6V3NiqoeYXJ2IyMX5FGSys7PZsWMHADNmzKBVq1YYhlH0el5eHu+9956CjIif8HgM5m3I4NWUdA7lngQgqVENxvVLpHVsdXOLExEpAZ+CjN1uZ/jw4Rw9ehSAu++++6zXb7/99tKvTkRK3eqdWUxI3swvGbkAxEaFMKZ3Ar2vqqsB7UTE7/gUZMLCwli1ahUAN954I4sWLSrTokSk9O0+ks+kr9NI2XQIgPCgAB7r2pR7OzQmOFATO4qIfypxH5nzhZjs7GyioqIuuyARKV05BS5mLN3Ge6t243IbWC1wZ/uGjOweT82wILPLExG5LCUOMj///DNTpkzh0KFDeDzezoEul4vs7Gx++eWXUi9QRC6Ny+3hwzV7mb5kK0cLXAB0io9mbN8E4uuEm1ydiEjpKHGQGT9+PLGxsTRr1ox9+/Zx3XXX8f777zNq1KiyqE9ESsgwDJamZzJxYRo7D+cD0Kx2GGP7JtC5eW2TqxMRKV0lDjLbtm3j//2//8f+/fuZOHEigwcPpk2bNowfP57BgweXRY0i4qO0X3OZmJzGiu1HAKhZzc7IHvHc0S6WAJvmRBKRyqfEQSYiIoLg4GBiY2PZtm0bAK1btyYjI6PUixMR32TmOZi6eCufrt+HxwC7zcrg6xszrEtTIoIDzS5PRKTMlDjIxMXF8dFHH/GXv/yF0NBQ0tLSsNvtemxTxAQOl5vZK3Yxc9l28p1uAPpeXY9nbmxBbFSoydWJiJS9EgeZESNG8Oijj3LdddfxwAMPcNttt2Gz2fjLX/5SFvWJyDkYhsH8nw4wZdEWMo6dAKBVg0ie65dIUmM9PSgiVUeJg0zbtm35/vvvCQwM5PbbbychIYG8vDyuu+66sqhPRM7ww56jTEjezIa9xwCoFxnM0ze24OZWMVg1saOIVDGXNBtcdnY2GRkZRdMU2O121q1bR7t27Uq1OBH53b7sAl5ZlM6Cn38FINRu49FOTXiwYxwhdg1oJyJVU4mDzFtvvcU//vGPs9ZbLBbS0tJKpSgR+V2ew8XM5TuYvWIXzkIPFgvcdk0so3rGUzsi2OzyRERMVeIgM3fuXN588026du2qDr4iZcjtMfhk3T6mfrOFI8edAFwbV5Nx/RK4MibS5OpERCqGEgeZgIAAOnfurBAjUob+u+0wE5PTSD+YB0BcrWo82yeBbgm19W9PROQ0JQ4yd911F9OmTeORRx4hLCysLGoSqbK2Z+YxMTmNZVsOAxAZEsiIbs24+4+NsAdoQDsRkTNd0jgyo0aNYvbs2We9pj4yIpcmO9/J9CVb+WDNXtwegwCrhXuubczwbk2pHmo3uzwRkQqrxEFm8uTJ3H///XTo0AGbTU9KiFyOk4Vu3l+5h9eXbiPPUQhAj8Q6jOndgrhotXiKiFxMiYNMXl6eJogUuUyGYZCy6SCTvk5nT1YBAIn1IhjXN4EOTWuZXJ2IiP8ocZDp0aMH33zzDT169CiLekQqvdT9OfwteTNrd2UDEB0exOhezRnQtgE2DWgnIlIiJQ4yDoeDESNG0KRJE6pXr17sCYr333+/VIsTqUwO5jiYkpLOvB+9E6wGBVgZckMcQzo1oVrQJY1NKSJS5ZX4t2fTpk1p2rRpWdQiUikVOAuZ9d1OZn2/A4fLA0D/NvUZ3as5MdVDTK5ORMS/lTjIPPbYY2VRh0il4/EYzNuQwasp6RzKPQlAUqMajOuXSOvY6uYWJyJSSfgcZF588UVefPFFxowZc95tJk2aVCpFifi71TuzmJC8mV8ycgGIjQphTO8Eel9VVwPaiYiUIp+DzKkJIkXk/HYfyWfS12mkbDoEQHhQAI91bcq9HRoTHKjhCkRESpvPQeall14C4I477qBVq1Znvf7999+XXlUifianwMWMpdt4b9VuXG4DqwXubN+Qkd3jqRkWZHZ5IiKVVon7yAwePJgff/yx2Lrjx48zYsQINmzYUGqFifgDl9vDh2v2Mn3JVo4WuADoFB/N2L4JxNcJN7k6EZHKz6cgs2fPHvr27Yvb7cYwDBISEs7apm3btqVenEhFZRgGS9MzmbgwjZ2H8wFoVjuMsX0T6Ny8tsnViYhUHT4FmUaNGvHZZ5+Rm5vLww8/zNtvv13s9aCgIOLj48ukQJGKJu3XXCYmp7Fi+xEAalazM7JHPHe0iyXApokdRUTKk8+3lk61wixYsIDY2Nii9cePH8dut2O3l3xiO7fbzX333Uf9+vWZPHnyBbdNSUlhypQpfPvttyU+j0hpyMxzMHXxVj5dvw+PAXablcHXN2ZYl6ZEBAeaXZ6ISJVU4v8+Op1Ohg0bBsA333zDH//4Rzp27MgPP/xQ4pO/8cYbrF+//oLbuFwu3n77bZ544gk9OSWmcLjcvLlsO11eXc7H67whpu/V9fh2VCfG9E5QiBERMVGJO/u+/PLL1K5dG8MwmDp1KsOHD6datWpMnjyZzz77zOfjrFq1isWLF9OzZ88Lbnf//fcTFBTEQw89xPz580tarsglMwyD+T8dYMqiLWQcOwFAqwaRPNcvkaTGUSZXJyIicAlBZsuWLfzzn/8kIyODvXv3cuedd1KtWjVee+01n4+RlZXF2LFjmTlzJnPnzr3gtq+++ip169Zl3rx5JS21GLfbfVn7n+94pX1cqRh+3HuUlxems2FfDgB1I4N5qmc8N11dD6vVor93ERHK9lro6zFLHGQKCwsxDIP//e9/XHnllYSFhZGdnU1QkG9jZXg8HkaPHs3gwYNp0aLFRbevW7duSUs8p9TU1FI5TnkdV8yRmV/I/0s9zv/2OQAItlno36IaN8VXI4hMfv450+QKRUQqHjOvhSUOMh06dOCvf/0r6enpPPDAA+zbt4+nnnqKzp07+7T/rFmzsNvtDBo0qKSnviwtW7bEZiu9kVXdbjepqamlflwxR56jkH9+t4M5K/fgLPRgscCf2zZgZPem1I4INrs8EZEKqSyvhaeOfTElDjJ/+9vfmDNnDklJSfTr14/MzEyuvPJKRo0a5dP+X375JZmZmSQlJQHgcHj/57tkyZKLdvy9HDabrUwCR1kdV8pHodvDp+v3M/WbLRw57gTg2riajOuXwJUxkSZXJyLiH8y8FpY4yAQFBeFyuZg7dy5vvvkm8+fPZ8OGDRw/fpyQkJCL7r9o0aJi3z/zzDMAF338WqS0/XfbYSYsSGPLoTwArqhVjWf7JNA9obYmdhQR8RMlfvx6xowZrFmzhtdff53AwEBq1apF3bp1mTBhQqkU1KZNGz2dJGVqe2Yeg99dy6DZa9lyKI/IkECe75dIyuM30COxjkKMiIgfKXGLzFdffcVHH31EnTreX/ihoaFMmjSJHj16XFIBZ7bEnG++pltvvZVbb731ks4hApCd72T6kq18sGYvbo9BgNXCPdc2Zni3plQPLfmAjiIiYr4SB5mCggKiorxjaJwaoC44OBirVUOzS8V0stDN+yv38PrSbeQ5CgHokViHMb1bEBcdZnJ1IiJyOUocZFq3bs0bb7zByJEji5rg//3vf9OyZctSL07kchiGQcqmg0z6Op09WQUAJNaLYFzfBDo0rWVydSIiUhpKHGTGjh3Lvffey3/+8x/y8/Pp06cP+fn5vPvuu2VRn8glSd2fw9+SN7N2VzYA0eFBjO7VnAFtG2Czqg+MiEhlUeIgExsbS3JyMsuXLycjI4O6devSuXNnwsLURC/mO5jjYEpKOvN+zAAgKMDKkBviGNKpCdWCSvzjLiIiFdwl/WYPCQmhd+/epV2LyCUrcBYy67udzPp+Bw6XB4D+beozuldzYqpffFgAERHxT/ovqvg1j8dg3oYMXk1J51DuSQCSGtVgXL9EWsdWN7c4EREpcwoy4rdW78xiQvJmfsnIBSA2KoQxvRPofVVdjQUjIlJFKMiI39l9JJ9JX6eRsukQAOFBATzWtSn3dmhMcKCmixARqUoUZMRv5BS4eH3pNt5ftRuX28BqgTvbN2Rk93hqhvk2+7qIiFQuCjJS4bncHj5cs5fpS7ZytMAFQKf4aMb2TSC+TrjJ1YmIiJkUZKTCMgyDpemZTFyYxs7D+QA0qx3G2L4JdG5e2+TqRESkIlCQkQop7ddcJiansWL7EQBqVrMzskc8d7SLJcCm6TBERMRLQUYqlMw8B1MXb+XT9fvwGGC3WRl8fWOGdWlKRHCg2eWJiEgFoyAjFYLD5Wb2il3MXLadfKcbgL4t6/FM7xbERoWaXJ2IiFRUCjJiKsMwmP/TAaYs2kLGsRMAtGoQyXP9EklqHGVydSIiUtEpyIhpfthzlAnJm9mw9xgA9SKDefrGFtzcKgarJnYUEREfKMhIuduXXcAri9JZ8POvAITabTzaqQkPdowjxK4B7URExHcKMlJu8hwuZi7fwewVu3AWerBY4LZrYhnVM57aEcFmlyciIn5IQUbKXKHbw6fr9zP1my0cOe4E4Nq4mozrl8CVMZEmVyciIv5MQUbK1H+3HWbCgjS2HMoD4Ipa1Xi2TwLdE2prYkcREblsCjJSJrZn5jExOY1lWw4DEBkSyIhuzbj7j42wB2hAOxERKR0KMlKqsvOdTF+ylQ/W7MXtMQiwWrjn2sYM79aU6qF2s8sTEZFKRkFGSsXJQjfvr9zD60u3kecoBKBHYh3G9G5BXHSYydWJiEhlpSAjl8UwDFI2HWTS1+nsySoAILFeBOP6JtChaS2TqxMRkcpOQUYu2c/7jzFhQRprd2cDEB0exOiezRlwTQNsGtBORETKgYKMlNivOSd4NWUL837MACAowMqQG+IY0qkJ1YL0IyUiIuVHVx3xWYGzkFnf7WTW9ztwuDwA9G9Tn9G9mhNTPcTk6kREpCpSkJGL8ngM5m3I4NWUdA7lngQgqVENxvVLpHVsdXOLExGRKk1BRi5o9c4sJiRv5peMXABio0IY0zuB3lfV1YB2IiJiOgUZOafdR/KZ9HUaKZsOARAeFMBjXZtyb4fGBAdqYkcREakYFGSkmJwCF68v3cb7q3bjchtYLXBn+4aM7B5PzbAgs8sTEREpRkFGAHC5PXy4Zi/Tl2zlaIELgE7x0Yztm0B8nXCTqxMRETk3BZkqzjAMlqZnMnFhGjsP5wPQrHYYY/sm0Ll5bZOrExERuTAFmSos7ddcJiansWL7EQBqVrMzskc8d7SLJcCmiR1FRKTiU5CpgjLzHExdvJVP1+/DY4DdZmXw9Y0Z1qUpEcGBZpcnIiLiMwWZKsThcjN7xS5mLttOvtMNQN+W9Ximdwtio0JNrk5ERKTkFGSqAMMwmP/TAaYs2kLGsRMAtGoQyXP9EklqHGVydSIiIpdOQaaS+2HPUSYkb2bD3mMA1IsM5ukbW3BzqxismthRRET8nIJMJbUvu4BXFqWz4OdfAQi123i0UxMe7BhHiF0D2omISOWgIFPJ5DlczFy+g9krduEs9GCxwG3XxDKqZzy1I4LNLk9ERKRUKchUEoVuD5+u38/Ub7Zw5LgTgGvjajKuXwJXxkSaXJ2IiEjZUJCpBP677TATFqSx5VAeAFfUqsazfRLonlBbEzuKiEilpiDjx7Zn5jExOY1lWw4DEBkSyIhuzbj7j42wB2hAOxERqfwUZPxQdr6T6Uu28sGavbg9BgFWC/dc25jh3ZpSPdRudnkiIiLlRkHGj5wsdPP+yj28vnQbeY5CAHok1mFM7xbERYeZXJ2IiEj5U5DxA4ZhkLLpIJO+TmdPVgEAifUiGNc3gQ5Na5lcnYiIiHkUZCq4n/cfY8KCNNbuzgYgOjyI0T2bM+CaBtg0oJ2IiFRxCjIV1K85J3g1ZQvzfswAICjAypAb4hjSqQnVgvTXJiIiAgoyFU6Bs5BZ3+1k1vc7cLg8APRvU5/RvZoTUz3E5OpEREQqFgWZCsLjMZi3IYNXU9I5lHsSgKRGNRjXL5HWsdXNLU5ERKSCUpCpAFbvzGJC8mZ+ycgFIDYqhDG9E+h9VV0NaCciInIBCjIm2n0kn0lfp5Gy6RAA4UEBPNa1Kfd2aExwoCZ2FBERuRgFGRPkFLh4fek23l+1G5fbwGqBO9s35PHu8dQKCzK7PBEREb9hapBxu93cd9991K9fn8mTJ59zm++++46///3v7Nu3j3r16vHUU0/RpUuXcq60dLjcHj5YvYfp327jWIELgE7x0Yztm0B8nXCTqxMREfE/pgaZN954g/Xr11O/fv1zvr57927++te/MnXqVDp37szixYt5/PHHWbx4MXXq1Cnnai+dYRgsTc9k4sI0dh7OB6BZ7TDG9k2gc/PaJlcnIiLiv0wLMqtWrWLx4sX07NnzvNv85z//ISkpie7duwPQp08f5s2bxyeffMLw4cNLdD63231Z9Z7veBc7bvrBPF5emM7/dmQBEFXNzshuTbktqQEBNmup1yUiIlJefL0WXs6xL8aUIJOVlcXYsWOZOXMmc+fOPe9227dvJz4+vti6pk2bkp6eXuJzpqamlnifyznuUYebj385ztJdJ/AAAVbo16watyZUo1pgNr+kZpdJPSIiIuWtrK6xvij3IOPxeBg9ejSDBw+mRYsWF9w2Pz+fkJDig8AFBwdTUFBQ4vO2bNkSm630ngRyu92kpqaedVyHy82c/+3mn9/tJN/pTZN9rqrLU73iiY0KLbXzi4iImO1818LSPPbFlHuQmTVrFna7nUGDBl1025CQEBwOR7F1DoeDatWqlfi8Nput1D/k049rGAbzfzrAlEVbyDh2AoBWDSJ5rl8iSY2jSv28IiIiFUVZXWN9Ue5B5ssvvyQzM5OkpCSAoqCyZMkS1q9fX2zb+Ph4Nm3aVGzd9u3bueqqq8qnWB/9sOcoE5I3s2HvMQDqRQbz9I0tuLlVDFZN7CgiIlJmyj3ILFq0qNj3zzzzDMA5H7+++eabeffdd1m4cCE9e/Zk8eLFrF27lrFjx5ZLrReTmV/I8I83kpx6EIBQu41HOzXhwY5xhNg1oJ2IiEhZs5pdwJnatGnD/PnzAWjSpAlvvvkms2bNol27dsycOZMZM2ZwxRVXmFwlfLZ+P8MXHSE59SAWC9yeFMvyJzvz127NFGJERETKiekj+57ZErNhw4Zi33fs2JGOHTuWZ0k+SU79FZcH/hgXxXP9ErkyJtLskkRERKoc04OMv5p2eyuWrfmJ/p2TCAjQxygiImKGCndryV/UCLUTVyNQs1OLiIiYSEFGRERE/JaCjIiIiPgtBRkRERHxWwoyIiIi4rcUZERERMRvKciIiIiI31KQEREREb+lICMiIiJ+S0FGRERE/JaCjIiIiPgtBRkRERHxWwoyIiIi4rcUZERERMRvKciIiIiI31KQEREREb+lICMiIiJ+S0FGRERE/JaCjIiIiPgtBRkRERHxWwoyIiIi4rcUZERERMRvKciIiIiI31KQEREREb+lICMiIiJ+S0FGRERE/JaCjIiIiPgtBRkRERHxWwoyIiIi4rcUZERERMRvKciIiIiI31KQEREREb+lICMiIiJ+S0FGRERE/JaCjIiIiPgtBRkRERHxWwoyIiIi4rcUZERERMRvKciIiIiI31KQEREREb+lICMiIiJ+S0FGRERE/JaCjIiIiPgtBRkRERHxWwFmF+C3tnzNFT+8jWV/Q6hWC0JrQmiU98+QqN+/Do4Ei8XsakVERColBZlLZF07i6gD38OBi2xosZ074BRbV/O0dVEQFAlWNZaJiIhcjILMJfLcMpOMpW/TICoEq+MYFGT9tmT/tmSBKx8MN+Qf9i6+stggpEbxwBMadY7Q89vXITUguLrCj4iIVDkKMpcqIobDjW+hfuvWYLOdexuXA05knxZwsn77/ox1p693HveGn4Ij3sVXFutp4cfH1h+FHxER8XMKMmUpMBgCYyAixvd9Ck+eEXpOtfQcPce6U+EnDwzP7+t9dSr8hEQVv7V1zltev60LqQ7W8wQ3ERGRcqYgU9EEBEFEPe/iq8KTcOLoGQHnQq0/R+FkbvHwk7XNx5NZfmv5ifK99UfhR0REyoiCTGUQEAThdb2LrwqdxcPPRVt/jsLJHMDwrj+RDVnbfTyZxRtmigWci7X+1FD4ERGRizIlyKxatYqpU6eyY8cOQkJCuPHGGxk9ejTBwcFnbTtv3jz+9a9/cejQIeLj43nyySdp166dCVVXMgF2CK/jXXx1Kvyc2e/nVCvPmS1CBdmnhZ+j3qUkgqufu2PzWetOCz82ZXMRkaqk3H/rZ2dnM2TIEF588UX+9Kc/ceTIER544AH+9a9/MXz48GLbfvvtt7zwwgu8/vrr3HDDDXz77bc89NBDzJs3j7i4uPIuXS4l/Lhdv4WcMzs2n/GE1+mtP44c776OY94le4fv5wuOPDvgXKj1J6QG2AJL8imIiEgFUu5BJioqipUrVxIWFoZhGBw7doyTJ08SFRV11rYLFiygX79+dOnSBYCePXvy6aef8n//93+MHj26vEuXS2ELhLDa3sVX7sLfW3jO6tx8ntYfxzHvvo4c75K90/fzBUWe3cfngq0/UQo/IiIVhCnt8GFhYQB06tSJQ4cOkZSUxK233nrWdm63m9DQ0GLrrFYrO3eW4CJ12rFK06njlfZxBbx9an4bN8dXnkI4caxY0LGc+K1vzwlv6LEUff3bnyeOYcHw3v46mQNHd/l8OiMootjYPsapgHPq+6KvT2sJUvgRkUqmLK+Fvh7TYhiGUepn95HD4SAnJ4cnn3ySoKAg3nnnnWKvJycn8/zzz/PWW2/Rtm1bli9fzsiRI2nXrh1z5szx6Rxut5uNGzeWQfXi9ww3NudxApw53sWV+9vXub+vO+NrmyvPG34ugTugGoX2iN+WSAoDf/vTHvn7Ovvv69z2CAyrwo+IVG2tW7fGdr7x2jD5qaXg4GCCg4MZPXo0f/7zn8nJySEyMrLo9b59+5Kdnc1zzz1HTk4OnTp1ol+/fpw4caLE52rZsuUFP4iScrvdpKamlvpxpWLzeNze21gF2UWtO96Wn+yidZaC7N8fff/tTwsGtsJ8bIX5BBX86vP5DHtYsVtaxhmtPMaZ01+ERHmfYhMRKQdleS08deyLKfcg8+OPP/Lss88yf/587HY7AE6nk8DAQEJCQopte/jwYTp27MigQYOK1t1222307NmzxOe12WxlEjjK6rhSQdlsEFgbwkvQ58fj9vbbOW+H56zTQs9pT4EZHizO497Rno/tAcCn6UftYReY0iLqjE7Qv4WgwLOfGBQR8ZWZ18JyDzLNmzfH4XDw2muvMWrUKA4fPswrr7zCwIEDi4LNKevWrWPSpEl8/PHH1KpVi48++ohdu3bRv3//8i5b5NJZT00cGgU09W0fj8fb8nPOx9rPM9jhiWzvIIdF4Wev7zUGVrvAE16nD3h42nqFHxGpAMo9yFSrVo133nmHl19+meuuu47w8HBuuukmhg0bBkCbNm146aWXuPnmm+nTpw87d+7k9ttvp6CggCuvvJL33nuPmjVrlnfZIuXLav09QNRs4ts+Ho+30/Lpj7VftPUn2zu3lysfcvIhp6ThJ6pkrT+BIRc/rohICZja2bc8nOrse7HOQhXluCLlyuPxTldxroBzodYf4xKfUAgMPfeUFsWmujhjnT304scVEVOU5bXQ12NrGFSRqsxq9U4fEVLd95Yfwygefnwd7NBTCK4C75K73/caA0J+CzY1zj+lxZmtP4GhYPGpR5GI+DkFGREpGYvFO4JycCRE+TjCdlH4yeasjs0XmurC44LCE97gU6LwE+zDhKZntP4o/Ij4JQUZESl7xcLPFb7tYxhwMu/crTwXav1xO6HQAbkZ3sVXtiAfJjQ9Y729msKPiMkUZESkYrJYIDjCu9Ro7Ns+huF9Yuv0gONL64/bCe6TkHfAu/jKZj/HlBYXaf2xhyn8iJQiBRkRqTwsFggK9y4lCj/552/lKbb+9PBz0huA8n71Lr6y2X2b0PT0fj9B4Qo/IuehICMiVZvFAkFh3qVGI9/2MQxvp+ViLT3na/05bX2hwxt+jh/0Lr6yBp5nQtMLtP4ERSj8SJWgICMiUlIWi7d/jL0aVG/o+37OgjNaec73uPup1p8j3vDjccHxQ97FV9aAMwLOaeHnfH1/FH7EDynIiIiUF3uod6ke6/s+zoJztPRcpPXHVeB93P1Sws+ZT3qd63H307cJjlT4EVMpyIiIVGSnwk9kA9/3KQo/2T60/vy2uPK94Sc/07v4ymI7R8A5x5QWxW57RXrHMBIpBQoyIiKVzaWEH9eJ87T0nGddQZY3/BhuyD/sXXxlsZ3Wz8fH1p/g6go/ck4KMiIi4p0HK7K+d/GVy3FGa88FJjQ99b3zuDf8FBzxLr6yWIt3cval9Ufhp0pQkBERkUsTGAyBMRAR4/s+hSfPEXBOe7T9XK0/zjzvzO6n1vvKYvWGmQsNanhmEAqp7p2xXvyGgoyIiJSfgCCIqOddfFV48rQBDH1p/cn2TolheLzfn8iGrG0+nszyW8uPDxOaFn1dXeHHRAoyIiJSsQUEQXhd7+KrQmfx0Zsv2vpzFE7mAMZp4We7jyez/Db56kWmtDi9BSi4Oth0CS4N+hRFRKTyCbBDeB3v4qtT4edE9jlaf86Y0PTUn0Xh56h3yd7h+/mKbnuda7DDc7T+hNRQ+DkHfSIiIiJwaeHH7fot5Pgwoemp1x053n0dx7xLicJPpG8TmhY9CVYDbIEl+RT8joKMiIjIpbIFQlht7+Ird+HvLTy+tv44jnn3deT8FoR2+n6+oMjij7lftPUnyq/Cj4KMiIhIebIFQFi0d/GVu9AbZi44pcUZ604cAwzv7a+TOXB0l+/nC4r8fVqLYq08ZzwCH1wdS6GjhB9A6VKQERERqehsAVCtlnfxlcftDTPnvOV1xoSmRd8fpXj42X3hsoCWQVHQKhVsYZfxBi+dgoyIiEhlZLVBtZrexVcet/fW1Vkdm8/d+mOcyOZkUF1CTHz8XEFGREREvKy23zsP0+yim3vcbrZs3Ehrm73sazsPjd0sIiIifktBRkRERPyWgoyIiIj4LQUZERER8VsKMiIiIuK3FGRERETEbynIiIiIiN9SkBERERG/pSAjIiIifktBRkRERPyWgoyIiIj4LQUZERER8VsKMiIiIuK3FGRERETEbwWYXUBZMwwDALfbXarHPXW80j6uiIiIvyjLa+GpY566jp+PxbjYFn7O6XSSmppqdhkiIiJyCVq2bIndbj/v65U+yHg8HgoLC7FarVgsFrPLERERER8YhoHH4yEgIACr9fw9YSp9kBEREZHKS519RURExG8pyIiIiIjfUpARERERv6UgIyIiIn5LQUZERET8loKMiIiI+C0FGREREfFbCjIiIiLitxRkRERExG8pyJzD888/z/33319s3fjx43nqqafYu3cvjzzyCO3bt6dLly5MmzYNp9MJwPHjxxk5ciTt27fnuuuu44EHHmDHjh1mvAUREZFSs3//fpo3b85nn31G165dueaaaxg8eDAHDx4EYMmSJdx66620bduWXr16MXfuXDweT7nUpiBzDgMHDmTVqlUcOnQI8E48mZycTJ8+fbjvvvto1qwZ33//PR9++CErV65kxowZAMyZM4fjx4/z3XffsWzZMqKjo/n73/9u5lsREREpNcuXL+eLL74gJSWFI0eOMHPmTFavXs3jjz/Ogw8+yNq1a5k6dSrvvvsu77//frnUpCBzDldffTVNmjRhwYIFgPcvLiwsjIKCApxOJ0888QRBQUHUq1ePESNG8MEHHwAQHBxMeno6X3zxBYcOHeLll1/mrbfeMvOtiIiIlJqHHnqIiIgIatWqRdeuXdm9ezfz5s2jW7du9OnTh4CAAK688koefvhhPv7443KpKaBczuKHbr31Vr744gseeOAB5s2bR//+/cnIyCA7O5t27doVbWcYBi6Xi6ysLB566CHsdjuff/4548ePJzY2llGjRtGzZ08T34mIiEjpqFWrVtHXAQEBGIZBVlYWCQkJxbZr0KABGRkZ5VKTgsx53HLLLUydOpUNGzbwv//9j+eff54ffviBhg0bsmjRoqLtjh8/TlZWFlFRUWzZsoWuXbty3333kZeXx4cffsjIkSNZvXo14eHhJr4bERGRslG/fn327t1bbN2+ffuIjo4ul/Pr1tJ51KxZk06dOjF+/HiSkpKIiYmhS5cu5Ofn88477+B0OsnNzeXpp59m5MiRWCwWPvvsM5566imysrIICwsjLCyM0NBQ7Ha72W9HRESkTAwYMIClS5fy9ddf43a72bx5M2+//TYDBgwol/MryFzArbfeyubNm4v+MsLCwpg7dy5r1qzhhhtuoHv37lit1qJ+ME888QSNGjWib9++tG3blnnz5jFz5kyCgoLMfBsiIiJlplWrVvzjH//g7bffJikpiccee4y//OUvPPLII+VyfothGEa5nMkPpaenM2jQIFasWKEwIiIiUgGpj8w5HD9+nAMHDjB9+nRuvfVWhRgREZEKSreWzuHgwYPcfvvt5OTkMHToULPLERERkfPQrSURERHxW2qREREREb+lICMiIiJ+S0FGRERE/JaCjIiIiPgtBRkRKXcnT57k4MGDZpchIpWAgoyIlLs777yTlStXml2GiFQCCjIiUu6OHj1qdgkiUkkoyIhIubr//vs5cOAAL7zwAuPHj2fTpk0MGjSIdu3a0bNnT+bOncup4a1mzJjB8OHDefLJJ0lKSuKGG27gtddeKzrWoEGDmDFjRtH3+/fvp3nz5uzfvx+A5s2bM2HCBNq3b18078vKlSsZOHAgSUlJ9O3bl/nz55fjuxeR0qYpCkSkXM2ZM4euXbvy2GOPcd1119G3b19GjhzJnDlz2LNnD0OHDiU4OJg77rgDgMWLFzN58mReeeUVVqxYwZAhQ+jWrRutW7f26Xx79+5l+fLluFwu0tPTefTRR3n11Vfp1q0bP/30E0OHDqVGjRp07NixDN+1iJQVtciIiGnmz59PkyZNuOuuuwgMDKRp06Y88MADfPDBB0XbNG7cmD/96U/YbDY6depEdHQ0u3fv9vkc/fr1IyQkhIiICD7++GO6detGz549sdlstG3blttuu63Y+UTEv6hFRkRMk5GRwaZNm0hKSipa5/F4sNlsRd9HR0cX2ycwMBCPx+PzOWrXrl3sfKtXry52PrfbTcOGDS+lfBGpABRkRMQ0devWpX379syePbto3dGjR8nPz/dpf6vVisvlKrbvmSwWS7Hz9e/fn/Hjxxety8zMRFPOifgv3VoSkXJnt9vJy8vjpptuYuPGjcyfP5/CwkIyMzN55JFHmDx5sk/HadKkCf/973/Jzc0lLy+Pt99++4LbDxw4kAULFrBixQo8Hg+7d+/m7rvvZs6cOaXxtkTEBAoyIlLuBg4cyLRp05g2bRrvvPMOn3zyCR06dOCWW24hLi7O5yAzZMgQatasSbdu3bjlllvo2rXrBbdv1aoVU6dOZerUqbRr1467776brl27MmrUqNJ4WyJiAouhNlURERHxU2qREREREb+lICMiIiJ+S0FGRERE/JaCjIiIiPgtBRkRERHxWwoyIiIi4rcUZERERMRvKciIiIiI31KQEREREb+lICMiIiJ+S0FGRERE/Nb/B8c9qM66YNjPAAAAAElFTkSuQmCC",
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Your answer here\n",
"# plot\n",
"sns.lineplot(data=anova_predictions,\n",
" y='estimate', x='tenure',\n",
" hue='gender')"
]
},
{
"cell_type": "markdown",
"id": "058ddc5f-eac2-4bc8-88fe-d3593a287224",
"metadata": {},
"source": [
"The ANOVA suggested we had the following results:\n",
"1. A main effect of gender (differences between men and women, ignoring tenure status)\n",
"2. A main effect of tenure (differences between tenured and non-tenured, ignoring gender)\n",
"3. An interaction, indicating the difference between one variable (e.g. gender) at one level of the other (say tenured) is different to the other (confusing!)\n",
"\n",
"Have the model make predictions, and to explore the main effects, use the `by` keyword to ignore one variable and the `hypothesis` keyword to check the differences."
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "748d0f0b-2dec-4771-ab4d-0ad0e44f1e7e",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"shape: (1, 8)
term
estimate
std_error
statistic
p_value
s_value
conf_low
conf_high
str
f64
f64
f64
f64
f64
f64
f64
"Row 1 - Row 2"
0.175352
0.060417
2.902376
0.003703
8.076918
0.056937
0.293766
"
],
"text/plain": [
"shape: (1, 8)\n",
"┌───────────────┬──────────┬───────────┬─────┬─────────┬──────┬────────┬───────┐\n",
"│ Term ┆ Estimate ┆ Std.Error ┆ z ┆ P(>|z|) ┆ S ┆ 2.5% ┆ 97.5% │\n",
"│ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- │\n",
"│ str ┆ str ┆ str ┆ str ┆ str ┆ str ┆ str ┆ str │\n",
"╞═══════════════╪══════════╪═══════════╪═════╪═════════╪══════╪════════╪═══════╡\n",
"│ Row 1 - Row 2 ┆ 0.175 ┆ 0.0604 ┆ 2.9 ┆ 0.0037 ┆ 8.08 ┆ 0.0569 ┆ 0.294 │\n",
"└───────────────┴──────────┴───────────┴─────┴─────────┴──────┴────────┴───────┘\n",
"\n",
"Columns: term, estimate, std_error, statistic, p_value, s_value, conf_low, conf_high"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Your answer here\n",
"# Main effects\n",
"# Gender\n",
"me.predictions(anova_model, newdata=anova_predmat, by='gender', hypothesis='pairwise')\n",
"\n",
"# Tenure\n",
"me.predictions(anova_model, newdata=anova_predmat, by='tenure', hypothesis='pairwise')"
]
},
{
"cell_type": "markdown",
"id": "39c2d04c-5d51-42eb-a11c-12bf50d0c9a6",
"metadata": {},
"source": [
"Now use the predictions to figure out the 'cause' of the interaction. There are a few ways to do this. You can compare men and women professors who are tenured, and see if that difference is significant, and then see if the difference between non-tenured professors is also significant. What do you observe?"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "50197fde-97b4-4ab9-826d-bdc0d4b5d3b2",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"shape: (2, 10)
gender
term
contrast
estimate
std_error
statistic
p_value
s_value
conf_low
conf_high
str
str
str
f64
f64
f64
f64
f64
f64
f64
"female"
"tenure"
"mean(yes) - mean(no)"
0.055172
0.08796
0.627241
0.530501
0.914573
-0.117227
0.227572
"male"
"tenure"
"mean(yes) - mean(no)"
-0.405876
0.082847
-4.899093
9.6280e-7
19.98626
-0.568254
-0.243499
"
],
"text/plain": [
"shape: (2, 10)\n",
"┌────────┬────────┬──────────────────────┬──────────┬───┬──────────┬───────┬────────┬────────┐\n",
"│ gender ┆ Term ┆ Contrast ┆ Estimate ┆ … ┆ P(>|z|) ┆ S ┆ 2.5% ┆ 97.5% │\n",
"│ --- ┆ --- ┆ --- ┆ --- ┆ ┆ --- ┆ --- ┆ --- ┆ --- │\n",
"│ str ┆ str ┆ str ┆ str ┆ ┆ str ┆ str ┆ str ┆ str │\n",
"╞════════╪════════╪══════════════════════╪══════════╪═══╪══════════╪═══════╪════════╪════════╡\n",
"│ female ┆ tenure ┆ mean(yes) - mean(no) ┆ 0.0552 ┆ … ┆ 0.531 ┆ 0.915 ┆ -0.117 ┆ 0.228 │\n",
"│ male ┆ tenure ┆ mean(yes) - mean(no) ┆ -0.406 ┆ … ┆ 9.63e-07 ┆ 20 ┆ -0.568 ┆ -0.243 │\n",
"└────────┴────────┴──────────────────────┴──────────┴───┴──────────┴───────┴────────┴────────┘\n",
"\n",
"Columns: gender, term, contrast, estimate, std_error, statistic, p_value, s_value, conf_low, conf_high"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Your answer here\n",
"me.predictions(anova_model, newdata=anova_predmat, hypothesis='b1=b2') # Tenured, NON significant\n",
"me.predictions(anova_model, newdata=anova_predmat, hypothesis='b3=b4') # Non-tenured, significant, males > females\n",
"\n",
"# Slopes\n",
"me.slopes(anova_model, newdata=anova_predmat, variables='tenure', by='gender')"
]
},
{
"cell_type": "markdown",
"id": "4a9d8c72-61ef-47d9-a4ad-a1439fd6bdce",
"metadata": {},
"source": [
"### e. ANCOVA done with marginal effects\n",
"Let us now add some complexity. ANCOVA is often described as an ANOVA 'adjusting' for another variable. We know it simply as a general linear model, with some kind of categorical predictor, and other continuous predictors that are also in the model. There can be as many categorical predictors and interactions between them as needed, as well as the continuous covariates.\n",
"\n",
"ANCOVA is a confusing and unnecessary term. Linear models are simpler, and here we will see how. \n",
"\n",
"First, carry out an ANCOVA with `pingouin` that looks at teaching evaluations between men and women (the categorical predictor), but adjusts for their beauty (the continuous covariate). Print the result. What does it tell you?\n"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "dd07d46a-4774-4d96-bc93-0f6348da1299",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
Source
\n",
"
SS
\n",
"
DF
\n",
"
F
\n",
"
p-unc
\n",
"
np2
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
gender
\n",
"
4.346745
\n",
"
1
\n",
"
15.055490
\n",
"
0.000120
\n",
"
0.031692
\n",
"
\n",
"
\n",
"
1
\n",
"
beauty
\n",
"
6.243877
\n",
"
1
\n",
"
21.626444
\n",
"
0.000004
\n",
"
0.044903
\n",
"
\n",
"
\n",
"
2
\n",
"
Residual
\n",
"
132.808865
\n",
"
460
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" Source SS DF F p-unc np2\n",
"0 gender 4.346745 1 15.055490 0.000120 0.031692\n",
"1 beauty 6.243877 1 21.626444 0.000004 0.044903\n",
"2 Residual 132.808865 460 NaN NaN NaN"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Your answer here\n",
"# ANCOVA in pingouin\n",
"pg.ancova(data=profs, dv='eval', between='gender', covar='beauty')"
]
},
{
"cell_type": "markdown",
"id": "f451aa14-65aa-44b0-bce4-b6cb7f7a1f66",
"metadata": {},
"source": [
"You should see that there are significant effects of both gender and beauty, but there's little information of use here. \n",
"\n",
"Fit a linear model that is equivalent to this ANCOVA, called `ancova_mod`. Print the summary."
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "2dac5a1f-0aeb-489e-aaf2-d59f67b7c76f",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"
OLS Regression Results
\n",
"
\n",
"
Dep. Variable:
eval
R-squared:
0.066
\n",
"
\n",
"
\n",
"
Model:
OLS
Adj. R-squared:
0.062
\n",
"
\n",
"
\n",
"
No. Observations:
463
F-statistic:
16.33
\n",
"
\n",
"
\n",
"
Covariance Type:
nonrobust
Prob (F-statistic):
1.41e-07
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
coef
std err
t
P>|t|
[0.025
0.975]
\n",
"
\n",
"
\n",
"
Intercept
3.8838
0.039
100.468
0.000
3.808
3.960
\n",
"
\n",
"
\n",
"
gender[T.male]
0.1978
0.051
3.880
0.000
0.098
0.298
\n",
"
\n",
"
\n",
"
beauty
0.1486
0.032
4.650
0.000
0.086
0.211
\n",
"
\n",
"
Notes: [1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
],
"text/latex": [
"\\begin{center}\n",
"\\begin{tabular}{lclc}\n",
"\\toprule\n",
"\\textbf{Dep. Variable:} & eval & \\textbf{ R-squared: } & 0.066 \\\\\n",
"\\textbf{Model:} & OLS & \\textbf{ Adj. R-squared: } & 0.062 \\\\\n",
"\\textbf{No. Observations:} & 463 & \\textbf{ F-statistic: } & 16.33 \\\\\n",
"\\textbf{Covariance Type:} & nonrobust & \\textbf{ Prob (F-statistic):} & 1.41e-07 \\\\\n",
"\\bottomrule\n",
"\\end{tabular}\n",
"\\begin{tabular}{lcccccc}\n",
" & \\textbf{coef} & \\textbf{std err} & \\textbf{t} & \\textbf{P$> |$t$|$} & \\textbf{[0.025} & \\textbf{0.975]} \\\\\n",
"\\midrule\n",
"\\textbf{Intercept} & 3.8838 & 0.039 & 100.468 & 0.000 & 3.808 & 3.960 \\\\\n",
"\\textbf{gender[T.male]} & 0.1978 & 0.051 & 3.880 & 0.000 & 0.098 & 0.298 \\\\\n",
"\\textbf{beauty} & 0.1486 & 0.032 & 4.650 & 0.000 & 0.086 & 0.211 \\\\\n",
"\\bottomrule\n",
"\\end{tabular}\n",
"%\\caption{OLS Regression Results}\n",
"\\end{center}\n",
"\n",
"Notes: \\newline\n",
" [1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
],
"text/plain": [
"\n",
"\"\"\"\n",
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: eval R-squared: 0.066\n",
"Model: OLS Adj. R-squared: 0.062\n",
"No. Observations: 463 F-statistic: 16.33\n",
"Covariance Type: nonrobust Prob (F-statistic): 1.41e-07\n",
"==================================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"----------------------------------------------------------------------------------\n",
"Intercept 3.8838 0.039 100.468 0.000 3.808 3.960\n",
"gender[T.male] 0.1978 0.051 3.880 0.000 0.098 0.298\n",
"beauty 0.1486 0.032 4.650 0.000 0.086 0.211\n",
"==================================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
"\"\"\""
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Your answer here\n",
"ancova_mod = smf.ols('eval ~ gender + beauty', data=profs).fit()\n",
"ancova_mod.summary(slim=True)"
]
},
{
"cell_type": "markdown",
"id": "42ec728f-3e19-4308-89f9-d1bed94dbcd1",
"metadata": {},
"source": [
"Now make a prediction about teaching evaluations for females and males. As the variable we want to control for is in the model (beauty), we don't need to make any predictions for it."
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "e5767d2b-27f9-407c-8596-2de46f959801",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"text/html": [
"
"
],
"text/plain": [
"shape: (1, 8)\n",
"┌───────────────┬──────────┬───────────┬──────┬─────────┬──────┬────────┬───────┐\n",
"│ Term ┆ Estimate ┆ Std.Error ┆ z ┆ P(>|z|) ┆ S ┆ 2.5% ┆ 97.5% │\n",
"│ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- ┆ --- │\n",
"│ str ┆ str ┆ str ┆ str ┆ str ┆ str ┆ str ┆ str │\n",
"╞═══════════════╪══════════╪═══════════╪══════╪═════════╪══════╪════════╪═══════╡\n",
"│ Row 1 - Row 2 ┆ 0.168 ┆ 0.0517 ┆ 3.25 ┆ 0.00115 ┆ 9.76 ┆ 0.0667 ┆ 0.269 │\n",
"└───────────────┴──────────┴───────────┴──────┴─────────┴──────┴────────┴───────┘\n",
"\n",
"Columns: term, estimate, std_error, statistic, p_value, s_value, conf_low, conf_high"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Your answer\n",
"# T-test style model\n",
"ancova_mod = smf.ols('eval ~ gender', data=profs).fit()\n",
"\n",
"# Predictions and contrast\n",
"me.predictions(ancova_mod, \n",
" hypothesis='pairwise',\n",
" newdata=me.datagrid(ancova_mod,\n",
" gender=['male', 'female'])\n",
" )"
]
},
{
"cell_type": "markdown",
"id": "53187ff3-129e-4454-846c-f3684e570af8",
"metadata": {},
"source": [
"### f. Knowledge of linear models gets you out of trouble\n",
"Following on from the last example, lets say we want to examine the interaction between gender and tenure status and control for beauty. Perhaps we wish to see whether our earlier ANOVA model stands up if we incorporate and control for beauty.\n",
"\n",
"First, try to fit one of these models in `pingouin`. Its another ANCOVA, but this time has two between factors."
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "6651ee2b-7edf-4569-a2b2-19395540440c",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": [
"hide-input"
]
},
"outputs": [],
"source": [
"# Your answer here\n",
"#pg.ancova(data=profs, dv='eval', between=['tenure', 'gender'], covar='beauty')"
]
},
{
"cell_type": "markdown",
"id": "fcc7e531-bd48-4731-9d7c-ecc44b052536",
"metadata": {},
"source": [
"If you did this correctly, you should see an error - the software doesn't support it!\n",
"\n",
"But we can easily fit a linear model to do this. Fit a model that has an interaction between gender and tenure, and has beauty as a predictor."
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "c9400fd1-5bec-4f6a-8f31-81f0b14de217",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"
OLS Regression Results
\n",
"
\n",
"
Dep. Variable:
eval
R-squared:
0.104
\n",
"
\n",
"
\n",
"
Model:
OLS
Adj. R-squared:
0.096
\n",
"
\n",
"
\n",
"
No. Observations:
463
F-statistic:
13.23
\n",
"
\n",
"
\n",
"
Covariance Type:
nonrobust
Prob (F-statistic):
3.32e-10
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
coef
std err
t
P>|t|
[0.025
0.975]
\n",
"
\n",
"
\n",
"
Intercept
3.8804
0.075
51.883
0.000
3.733
4.027
\n",
"
\n",
"
\n",
"
gender[T.male]
0.4890
0.105
4.650
0.000
0.282
0.696
\n",
"
\n",
"
\n",
"
tenure[T.yes]
0.0076
0.087
0.087
0.930
-0.164
0.179
\n",
"
\n",
"
\n",
"
gender[T.male]:tenure[T.yes]
-0.3668
0.121
-3.027
0.003
-0.605
-0.129
\n",
"
\n",
"
\n",
"
beauty
0.1289
0.032
4.032
0.000
0.066
0.192
\n",
"
\n",
"
Notes: [1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
],
"text/latex": [
"\\begin{center}\n",
"\\begin{tabular}{lclc}\n",
"\\toprule\n",
"\\textbf{Dep. Variable:} & eval & \\textbf{ R-squared: } & 0.104 \\\\\n",
"\\textbf{Model:} & OLS & \\textbf{ Adj. R-squared: } & 0.096 \\\\\n",
"\\textbf{No. Observations:} & 463 & \\textbf{ F-statistic: } & 13.23 \\\\\n",
"\\textbf{Covariance Type:} & nonrobust & \\textbf{ Prob (F-statistic):} & 3.32e-10 \\\\\n",
"\\bottomrule\n",
"\\end{tabular}\n",
"\\begin{tabular}{lcccccc}\n",
" & \\textbf{coef} & \\textbf{std err} & \\textbf{t} & \\textbf{P$> |$t$|$} & \\textbf{[0.025} & \\textbf{0.975]} \\\\\n",
"\\midrule\n",
"\\textbf{Intercept} & 3.8804 & 0.075 & 51.883 & 0.000 & 3.733 & 4.027 \\\\\n",
"\\textbf{gender[T.male]} & 0.4890 & 0.105 & 4.650 & 0.000 & 0.282 & 0.696 \\\\\n",
"\\textbf{tenure[T.yes]} & 0.0076 & 0.087 & 0.087 & 0.930 & -0.164 & 0.179 \\\\\n",
"\\textbf{gender[T.male]:tenure[T.yes]} & -0.3668 & 0.121 & -3.027 & 0.003 & -0.605 & -0.129 \\\\\n",
"\\textbf{beauty} & 0.1289 & 0.032 & 4.032 & 0.000 & 0.066 & 0.192 \\\\\n",
"\\bottomrule\n",
"\\end{tabular}\n",
"%\\caption{OLS Regression Results}\n",
"\\end{center}\n",
"\n",
"Notes: \\newline\n",
" [1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
],
"text/plain": [
"\n",
"\"\"\"\n",
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: eval R-squared: 0.104\n",
"Model: OLS Adj. R-squared: 0.096\n",
"No. Observations: 463 F-statistic: 13.23\n",
"Covariance Type: nonrobust Prob (F-statistic): 3.32e-10\n",
"================================================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"------------------------------------------------------------------------------------------------\n",
"Intercept 3.8804 0.075 51.883 0.000 3.733 4.027\n",
"gender[T.male] 0.4890 0.105 4.650 0.000 0.282 0.696\n",
"tenure[T.yes] 0.0076 0.087 0.087 0.930 -0.164 0.179\n",
"gender[T.male]:tenure[T.yes] -0.3668 0.121 -3.027 0.003 -0.605 -0.129\n",
"beauty 0.1289 0.032 4.032 0.000 0.066 0.192\n",
"================================================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
"\"\"\""
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Your answer here\n",
"# More complex ANCOVA\n",
"ancova2 = smf.ols('eval ~ beauty + gender * tenure', data=profs).fit()\n",
"ancova2.summary(slim=True)"
]
},
{
"cell_type": "markdown",
"id": "b201bade-ae5e-48ba-a369-87af55f8fd05",
"metadata": {},
"source": [
"Once you have this model, use it to make predictions about gender and tenure as before, and work out the interaction effects. Are they the same as before?"
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "af123534-3d80-42b2-b298-86ca46d7aff9",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"text/html": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Your answer here\n",
"sns.lineplot(data=me.predictions(ancova2, newdata=anova_predmat),\n",
" y='estimate', x='tenure',\n",
" hue='gender')"
]
},
{
"cell_type": "markdown",
"id": "96b774ec-7b88-48b5-abb9-4b92d245367c",
"metadata": {},
"source": [
"### g. Interpreting complex interactions with marginal effects\n",
"If you've completed the above exercises, you've mastered 99% of the statistics used in basic psychology, and learned how to do it from a much clearer perspective. Lets now build knowledge of how to interpret an even more complex model. \n",
"\n",
"Lets suppose that, rather than controlling for beauty's influence on teaching evaluations for tenured and non-tenured males and females, you want to know whether it influences evaluations at these combinations. That is, you might wish to see how less attractive males are evaluated before and after tenure, and whether this change is different for females, who are typically judged more harshly on their looks. \n",
"\n",
"To do this, you will need an interaction between gender, beauty, and tenure. Fit a model that does this, and call it `three_interact`. Print the summary."
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "964e9828-d394-4fba-9919-214f51e87a9b",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"
OLS Regression Results
\n",
"
\n",
"
Dep. Variable:
eval
R-squared:
0.110
\n",
"
\n",
"
\n",
"
Model:
OLS
Adj. R-squared:
0.097
\n",
"
\n",
"
\n",
"
No. Observations:
463
F-statistic:
8.055
\n",
"
\n",
"
\n",
"
Covariance Type:
nonrobust
Prob (F-statistic):
3.03e-09
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
coef
std err
t
P>|t|
[0.025
0.975]
\n",
"
\n",
"
\n",
"
Intercept
3.8601
0.076
51.031
0.000
3.711
4.009
\n",
"
\n",
"
\n",
"
gender[T.male]
0.5076
0.107
4.741
0.000
0.297
0.718
\n",
"
\n",
"
\n",
"
tenure[T.yes]
0.0275
0.088
0.312
0.755
-0.146
0.201
\n",
"
\n",
"
\n",
"
gender[T.male]:tenure[T.yes]
-0.3781
0.122
-3.100
0.002
-0.618
-0.138
\n",
"
\n",
"
\n",
"
beauty
0.0006
0.080
0.008
0.994
-0.156
0.157
\n",
"
\n",
"
\n",
"
gender[T.male]:beauty
0.1362
0.124
1.094
0.274
-0.108
0.381
\n",
"
\n",
"
\n",
"
beauty:tenure[T.yes]
0.1301
0.099
1.315
0.189
-0.064
0.325
\n",
"
\n",
"
\n",
"
gender[T.male]:beauty:tenure[T.yes]
-0.0934
0.146
-0.640
0.522
-0.380
0.193
\n",
"
\n",
"
Notes: [1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
],
"text/latex": [
"\\begin{center}\n",
"\\begin{tabular}{lclc}\n",
"\\toprule\n",
"\\textbf{Dep. Variable:} & eval & \\textbf{ R-squared: } & 0.110 \\\\\n",
"\\textbf{Model:} & OLS & \\textbf{ Adj. R-squared: } & 0.097 \\\\\n",
"\\textbf{No. Observations:} & 463 & \\textbf{ F-statistic: } & 8.055 \\\\\n",
"\\textbf{Covariance Type:} & nonrobust & \\textbf{ Prob (F-statistic):} & 3.03e-09 \\\\\n",
"\\bottomrule\n",
"\\end{tabular}\n",
"\\begin{tabular}{lcccccc}\n",
" & \\textbf{coef} & \\textbf{std err} & \\textbf{t} & \\textbf{P$> |$t$|$} & \\textbf{[0.025} & \\textbf{0.975]} \\\\\n",
"\\midrule\n",
"\\textbf{Intercept} & 3.8601 & 0.076 & 51.031 & 0.000 & 3.711 & 4.009 \\\\\n",
"\\textbf{gender[T.male]} & 0.5076 & 0.107 & 4.741 & 0.000 & 0.297 & 0.718 \\\\\n",
"\\textbf{tenure[T.yes]} & 0.0275 & 0.088 & 0.312 & 0.755 & -0.146 & 0.201 \\\\\n",
"\\textbf{gender[T.male]:tenure[T.yes]} & -0.3781 & 0.122 & -3.100 & 0.002 & -0.618 & -0.138 \\\\\n",
"\\textbf{beauty} & 0.0006 & 0.080 & 0.008 & 0.994 & -0.156 & 0.157 \\\\\n",
"\\textbf{gender[T.male]:beauty} & 0.1362 & 0.124 & 1.094 & 0.274 & -0.108 & 0.381 \\\\\n",
"\\textbf{beauty:tenure[T.yes]} & 0.1301 & 0.099 & 1.315 & 0.189 & -0.064 & 0.325 \\\\\n",
"\\textbf{gender[T.male]:beauty:tenure[T.yes]} & -0.0934 & 0.146 & -0.640 & 0.522 & -0.380 & 0.193 \\\\\n",
"\\bottomrule\n",
"\\end{tabular}\n",
"%\\caption{OLS Regression Results}\n",
"\\end{center}\n",
"\n",
"Notes: \\newline\n",
" [1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
],
"text/plain": [
"\n",
"\"\"\"\n",
" OLS Regression Results \n",
"==============================================================================\n",
"Dep. Variable: eval R-squared: 0.110\n",
"Model: OLS Adj. R-squared: 0.097\n",
"No. Observations: 463 F-statistic: 8.055\n",
"Covariance Type: nonrobust Prob (F-statistic): 3.03e-09\n",
"=======================================================================================================\n",
" coef std err t P>|t| [0.025 0.975]\n",
"-------------------------------------------------------------------------------------------------------\n",
"Intercept 3.8601 0.076 51.031 0.000 3.711 4.009\n",
"gender[T.male] 0.5076 0.107 4.741 0.000 0.297 0.718\n",
"tenure[T.yes] 0.0275 0.088 0.312 0.755 -0.146 0.201\n",
"gender[T.male]:tenure[T.yes] -0.3781 0.122 -3.100 0.002 -0.618 -0.138\n",
"beauty 0.0006 0.080 0.008 0.994 -0.156 0.157\n",
"gender[T.male]:beauty 0.1362 0.124 1.094 0.274 -0.108 0.381\n",
"beauty:tenure[T.yes] 0.1301 0.099 1.315 0.189 -0.064 0.325\n",
"gender[T.male]:beauty:tenure[T.yes] -0.0934 0.146 -0.640 0.522 -0.380 0.193\n",
"=======================================================================================================\n",
"\n",
"Notes:\n",
"[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
"\"\"\""
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Your answer here\n",
"# three variable interaction\n",
"three_interact = smf.ols('eval ~ gender * beauty * tenure', data=profs).fit()\n",
"three_interact.summary(slim=True)"
]
},
{
"cell_type": "markdown",
"id": "b07f0191-41fa-47d4-9c4b-788077a41835",
"metadata": {},
"source": [
"Confusion abounds looking at the coefficients. Lets make sense of this by asking for predictions from the model. Generate a data grid that asks for beauty to be evaluated at [-2, 0, 2] (that's -2 standard devs below average, average, and 2 above, as the variable is scaled so by the authors), for tenured and non-tenured males and females."
]
},
{
"cell_type": "code",
"execution_count": 22,
"id": "ba8978f1-5faf-43c9-a7e7-d8416d1c9481",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"text/html": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Your answer here\n",
"# Predictions\n",
"three_preds = me.predictions(three_interact, newdata=predmat)\n",
"\n",
"# Plot to show interaction pattern overall - lots of ways to do this, e.g.\n",
"sns.lineplot(data=three_preds,\n",
" x='beauty', \n",
" y='estimate',\n",
" style='gender',\n",
" hue='tenure')\n",
"\n",
"# Or\n",
"sns.relplot(data=three_preds,\n",
" x='beauty', y='estimate',\n",
" style='tenure', col='gender',\n",
" kind='line')\n"
]
},
{
"cell_type": "markdown",
"id": "9602f4ba-6845-414b-96f1-fb59f5e2a9e9",
"metadata": {},
"source": [
"If you have visualised it correctly, you should see the general pattern that male evaluations increase with beauty, but they are lower with tenure. Females on the other only show a positive beauty association *with* tenure, and no association without it.\n",
"\n",
"Now, are those differences meaningful? To test that we need to make a decision about how we want to evaluate our interaction. That depends on the question, since there are many ways to interpret interactions of this complexity.\n",
"\n",
"Do this in steps. First, is the association between beauty and evaluations different for females and males? "
]
},
{
"cell_type": "code",
"execution_count": 24,
"id": "855c8231-82f0-462e-b134-4cbf9e439ec3",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": [
"hide-input"
]
},
"outputs": [
{
"data": {
"text/html": [
"